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SUMMARY 


This i* the final report for NASA Contract NAS8-3313S. This report 
presents the results of Investigations on the various aspects of multiple 
scattering effects on visible and infrared laser beams transversing dense 
fog oil aerosols contained in a chamber (4* x 4* x 9'). In connection 
with these investigations the data for the time-history of aerosol 
optical depth, and the aerosol size distribution and mass concentration 
measurements were provided by Dr. Gerald C. Holst, Chemical Systems 
Laboratory. The report briefly describes: (i) the experimental details 
and measurements (Section 2) ; (il) analytical representation of the 
aerosol size distribution data by two analytical models— the regularized 
power law distribution and the Inverse modified gamma distribution 
(Section 3); (ill) retrieval of aerosol size distributions from multi- 
spectral optical depth measurements by two methods — the 2- and 3-parameter 
fast table search methods and the nonlinear least squares method (Section 4) : 
(iv) modeling of the effects of aerosol microphysical (coagulation and 
evaporation) and dynamical processes (gravitational settling) on the 
temporal behavior of aerosol size distribution, and hence on the extinc- 
tion of 4 laser beams with wavelengths - 0.44, 0.6328, 1.15 and 3.39 ym 
(Section 5) ; (v) and the exact and approximate formulations for four 
methods — Tam and Zardecici, Dolin-Fante, Arnush-Stotts , and our exact 
method — for computing the effects of multiple scattering on the trans- 
mittance of laser beams in dense aerosols, all of which are based on the 
solution of the radiative transfer equation under the small angle approxi- 
mation (Section 6) . In addition, the report identifies problem areas in 
which further research needs to be performed in the near future. 


a 



1 . INTRODUCTION 


Whttnavar a beam of radiation pasaea through an aeroaol mediiim, ita energy 
alwaya decreaaea during tranait. The loaa of energy ia due to either 
acattering or abaorption or both in the medium. The term extinction ia 
defined aa the aum of abaorption and acattering. The tranamittance T for a 
parallel, monochromatic beam of radiation of wavelength X after paaaage 
through an homogeneoua aeroaol layer of length L and extinction coefficient 

in the aingle acattering (SS) apprcximBtion, ia given by Bouguer'a 
-TtX) 

law, T - e , where the optical depth T(X) ■■3 (X) L. 

In SS, the acattered photon reaching the detector auf fera «. .tly one 
acattering event; whereaa, in multiple scattering (MS), it auffera more than 
one acattering event. Thua, the N order acattering refer a to a caae when a 
photon ia. jicattered and reacattered N timea before reaching the detector. 

The basia of the SS theory ia that i.f one knowa the scattering properties 
of an individual particle, then the scattering, effect of N similar particles 
is simply N times that of a single particle (Ref. 1 ) . Such a simple direct 
proportionality to N particles is what makes the SS theory so simple. 

In the case of MS, this proportionality to N particles no longer holds, 
thereby making che theory very complicated. Therefore, except for a few 
researchers (Refs. 1-R) , not much work has been done on experimental investiga- 
tions on MS; even though relatively more work has been done on the theoretical 
aspects of MS (Refs. 9-19). Therefore, a coordinated research program was 
undertaken between the Institute for Atmospheric Optics and Remote Sensing 
(IFAORS) , NASA-Marshall Space Flight Center (MSFC) and Army Aberdeen Chemical 
System Laboratory (CSL) to implement systematic research effort in the 
theoretical and experimental aspects of MS in aerosols in controlled laboratory 
environment. The theoretical research was performed at the Institute and the 
experimental work, at the NA! k-MSFC and CSL. Important features of such a 
theoretical and experimental approach to MS investigations was that the 
numeric^;! results based on theoretical MS models could be compared with 
experimental results and vice versa. 
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In this program, it was decided to start with the simple experiments on 
measurement of attenuation of visible and IR laser beams traversing through 
fog oil aerosols contained in a closed char.iber. Simultaneous measurements 
of optical depth for three or four wavelengths were made as a function of 
time. In addition, measurements of the aerosol slse distribution (SD) were made 
at two different times during such an experiment. The experimental details 
and the measurements are described in Section 2. 

In order to get a better understanding of the experimental SD 
data, the first step was to represent the data with some commonly used 
analytic model (s) (Ref. 21). The results are described in Section 3. 

The experimental results for T(A,t), provided us with data for three or 
four independent measurements of t, corresponding to the three or four wavelengths, 
times t during the experiments, which were conducted for time periods of 
about so to 150 min. From these three independent measurements, it was 
possible to retrieve three un)cnown parameters characterizing the aerosol 
size distribution at time t. There are several methods for retrieving aerosol 
size distribution from t versus X data (Ref. 22) . But, since the ultimate goal 
of our work was to develop a software that will provide on-line retrievals of 
the aerosol size distribution from laser beam extinction measuremet^ts, we 
worked on developing fast, approximate search methods, and compared the results 
with those obtained using a numerical inversion scheme, such as the Nonlinear 
Least Squares (NLLS) method. Validation of the results was conducted by 
checking them against the actual size distribution measurements. We developed 
the two-parameter and the three-parameter search methods, which will be dis- 
cussed in Sections 4.1, 4.2, and 4.3 of this report. An error analysis of 
the retrieved size distribution results is given in Section 4.4. 

The size distribution measurements and the retrievals show that as the 
time increases, the peak of the aerosol size distribution tends to shift toward 
larger radii and lower number densities. This means that processes other than 
simple sedimentation are also present. In order to explain the time behavior 
of the aerosol size distribution, we have to perform modeling tnd 3ensitivity 
studies, which involve developing models of aerosol microphysical and dynamical 
processes — considered separately or in combinations — to study how each process 
effects the time behavior of the size distribution and the optical depth 
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at tha wavalangths involvad. Tha rasulta of a paramatric modeling atudy to 
Invaatigate tha aeparata and combined af facta of aadimantation, coagulation 
and evaporation/growth are diacuaaad in Sactiona 5.1 and 5.2. Th«) theoretical 
atudy of multiple acattering affects of laser beams in danse aaroscls is 
discusaed in Section 6. 
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2 . EXPERIMENTAL DETAILS AND MEASUREMENTS 


Tha expariiwnta on laser beam propagation through dansa aerosols wera 
conducted at the Chemical Systems Laboratory (CSL) by Dr. Gerald Holst. 

Aerosc Is were placed Inside a closed chamber » which is essentially a 
wooden box (approximately 9* x S' x 4') blackened flat black on the 
inside and having entrance and exit apercures on opposite sides of the 
bo'A for the passage of laser beams (Fig. 1) . The experimental geometry 
is schematically illustrated in Fig. 2. 

The laser beams used in the experiments weret A « 441.6 nm (He-Cd) i 
612.J nm (He-Ne); 1.15 (Jm (He-Ne) ; and 3.30 (im (He-Ne) . The physical and 
optical properties of fog oil, used for aerosol generation, are given in 
Table i. 

Measurements of the time variation of optical depth for 3 or 4 
wavelengths simultaneously, and of size distributions were made in a 
series of experiments. Data for three sets of these experiments 
performed on November 30, 1978, January 21, 1979 and February 22, 1979 
were provided to us by Dr. G. Holst. Those are summarized in Table 2. 
Keasurements of x(A,t) and n(r) are described next. 

Multiwavelength Optical Depth vs. "*ime Curves ; Typically, with the 
box completely en^ty of fog oil aerosols, the laser beams were turned on, 
so that their beams after traversing through t)\e empty box impinged on 
lens-pinhole-type detector systems (Fig. 2b) , having very narrow fields of 
view, placed across the box opposite the lasers. The divergence angles of 
the laser beam was 'v 2 millirad and the acceptance angle of the detector 
field of view was 13 millirad. Narrow angles were necessary to ensure 
that the errors due to forward scattering were negligibly small. The 
intensities were recorded on chart recorders which were adjusted to full 
scale, representing 100 percent transmission, when the box was empty. 

The cheunber was then filled with fog oil aerosols by using the evaporation/ 
condensation technique, which involves letting a certain amount of fog oil 
vapor inside the box through an inlet, and stirring the air in the box by 
turning on three or four upward-directed fans located at the bottom of the 
box. The chart recorders, one for each of the laser beams, continued to 
record the direct beeun intensities, simultaneously. During the filling of the 
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FIGURE 1. Schematic illustration of the experimental apparatus. 
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(b) 


FIGURE 2, Schematic illustration of (a) the experimental geometry, 
and (b) the detector system. 
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TABLE 1. PROPEfVrZBS OP FOG OZL (SGP «2) 


Holal Avar age 

Boiling Point « 354*C 

Surfaca Tanaion <25*0 »• 34 dynes/cm 

Spacific Gravity (16*0 ■ 0.9218 gm/an^ 

Molecular Weight ■* 7'5 

Vapor PresBure (99*0 ■ 0.028 iwn of Hg 


A (pm) 

Refractive Index 

0.4416 

1.5077-i 

(2.339 X lO"®)* 

0.6328 

1.5077-i 

(9.527 X iO"®)« 

1.15 

1.52 - i 

(1.378 X lO"®)* 

3.39 

1.44 - i 

(0.101)** 

10.59 

1.48 - i 

(0.163)** 


^Obtained by Abbe Refractometer . Private communication by Dr. G. Holst 
**Ref. 1: Dr. M. R. Querry's Valuaa. Private communication. 





TABI^ 2. EXPERIMENTAL DATA SETS* FOR T VS. t and n(r) 


Meaauronent 

Dates 

Experimental Data 

Assumed 

Sat 



Taken 

Error 

Sat 1 

11/30/78 

(a) 

T vs. t for X - 0.6328* 1.15* and 
3.39 Dm 

± 3' 



(b) 

n(r) at time t - 0.5 and 51.25 min 




(c) 

mass concentration at several 

times during the experiment 


Set 2 

1/22/79 


T vs. t for X ■ 0.6328* 1.15* and 
3.39 Dm 


Set 3 

2/22/79 

(a) 

T VS. t for X « 0.4416* .6328* 1.15* 
and 3 . 39 Vim 

♦ 3% 



(b) 

n(r) at time t « 2.5 and 62.00 min 




(c) 

masc concentration M^, at several 
times during the experiment 


* 

Provided by Dr. 

G. Holst* 

CSL 
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chMibvrs with tog oil aeroBola* th« intanBltiss of th« direct beams continued 
to decrease, tfhen the aerosol appeared to be uniformly mixed, the fan was 
turned off, the time beln 9 referred to as t ■ 0 min. Thereafter, the Intensity 
of 0.4416 and 0.6328 Mm beams started to Increase, while that of the other two 
laser beams (X«l.l5 pm and 3.39 pro) first continued to decrease gradually and 
after a certain time started to Increase. From the chart-recorder Intensity 
measurements. It was easy to obtain t*..poral records of the optical depths, 1, 
for the laser beams by the use of Bouguer <Lambert-Beer) transmission law-, 

l(X.t) _ -T(X,f) 

uw 

where t is the time In minutes. The plots of ‘((X,t) vs. t curves for 
X - 0.6328, 1.15 and 3.39 pm for Set I (November 30, 1978) are shown In 
Fig. 3. They show that for X - 1.15 and 3.39 pm, T first Increased for nbout 
20 and 30 min, respectively, and then started to decrease, Indicating that 
the size distribution was changing with time. Similar plots of T(X,t) vs t 
for Sets 2 and 3 are shown in Figs. 4 and 5 respectively. 

In addition, experiments were performed with* two 0.6328 pm lasers placed 
at two different heights (see Fig. 2a) to ma)ce T vs. t measurements. In order 
to investigate if there was any difference In the T vs. t curves. Results 
for one such experiment are shown in Fig. 4 with curves mar)(ed upper and 
lower, which were close to within experimental error. 

Aerosol Size Distribution Measurements ! In addition, aerosol size 
distribution measurements were made with a ten-stage Andersen in^actor twice 
during the transmission experiments for Sets I and 3. For Set 1, the measure- 
ments were made one at the beginning (t^ ■ 0.5 min) and the second time near 

the end (t_ ■ 51. 25 min) and are shown in' Fig. 6. From dN/d log d data, we 

2 2 
derived the area artd volume size distributions by calculating d (dN/d log d) 

and d^(dN/d log d) . The latter are useful representations as they give us 

some idea of the scattering effects and the mass distribution of aerosols. 

The three types of representations of aerosol size distribution are plotted 

for Set 1 and times t •« 0.5 and 51.25 min in Figs. 7 wd 8 on log-log graphs, 

which are identical in scale to the graphs in Ref. 21. For Set 3, the size 

distribution (dN/d log d) measurements were made at t •= 2.5 min and 62.0 min, 

and the data is shown in Fig. 9. The data for mass concentration measured 

at various times during the Set 1 and Set 3 experiments are shown in Fig. 10. 
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FIGURE 4. Set 2 data for T vs. t for \ » 0.6328 yjn laser beam located 
at 30 cm and 127 cm levels from the top of the box, and for 
\ - 1.15 and 3.39 yra laser beams at 63.5 cm level (see Fig. 2a) 
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FIGURE 6. Set 1 data for the size distribution at times t » 0.5 and 
51.25 min. 
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AERODYNAMIC DIAMETER d, ym 

FIGURE 9. Set 3 data for size distributions at times t " 2.5 and 62 min. 
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3 . ANALYTIC REPRESENTATION OF AEROSOL 
SIZE DISTRIBUTION DATA 


Th« Andersen Impactor data for the aerosol sice distribution, dN/d 
measured at times t > 0.5 min and 51.25 min (Pig. 6) for Set 1 was converted 
to dN/dr, i.e. n(r), where r “ d/2, and plotted on a log-log graph identical 
in size to the graphs in the parametrized graphical catalog for size distribu- 
tion models (Ref. 21). By overlaying and matching the experimental data with 
the appropriate catalog graphs belonging to the selected analytic model, such 
as. Regularized Power Law Distribution (RPLD) , lnverr,« Modified Gamma 
Distribution (IMGD), or Log Normal Distribution (LND) , we obtained the first 
guess estimates for the model parameters p^, i « 1, 2, ... Then by inputing 
these values as 1st guess estimates of the model parameters into a non-linear 
least squares (NT1.S) program, we could quickly obtain the values of the p^ 
that provide the best fit to the experimental data. The n(r) data was thus 
fitted, in the least squares sense, with the following analytic models; 

RPLD, IMGD, LND, and two term models such as the IMGD plus power law and 
the LND plus power law distribution. From among these models, two gave the 
best fit to the experimental data. They are the Regularized Power Law 
Distribution (RPLD) 


n(r) 


P2 


(r/P2> 


P3-I 


P3 P4 

(1 (r/pj) ■*) ’ 


(la) 


and the Inverse Modified Gamma Distribution (IMGD) 


P4 

-P3/(r ^) 

n(r) - p, r (lb) 

^2 

r 

The results of the fits to these two models are shown ir Figs. 10-13 and 
discussed as follows: 

(i) Figure 11 (Data at t • 0.5 min, RPLD model).- The RPLD model gives 
good fitc to the data for r < 0.1 pm or r > 2.0 pm; but in the 0.1 - 2.0 pm 
range the fit is not as good. (Circles denote the experimental data; squares. 
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n(r) 




SIZE DISTRIBUTION 
SET 1 DATA 
TIME > 0.5 MIN 

MODEL 2 - RPL 


O MEASURED DATA 

O INITIAL ESTIMATES 
- 2.307E + 07 

P 2 B .250 

pj » 6.00 

P 4 - 1.67 
O NLLS BEST FIT 
Iterations “ 5 
p^ - 1.916E + 07 
P2 • .'^21 

Pj ■ 6.17 
P 4 « 1.75 


RADIUS r. pm 

FIGURE 11. Analytic representation of Set 1, t • 0.5 min, size 

distribution data using a Regularized Power Law (RPL) Model 



the initial estimatea of the model parameters p^. i - 1, 2, 3, and 4 | diamonds, 
the best fit values of p^, obtained after convergence is reached in 7 iterations.) 

(ii) Figure 12 (Date at t - 0.5 min, IMGD modal}: Good fits to data 

for r < 0.2 pm wd r > 0.7 pm. The fit around the peak is somewhat poor. 

(Hi) Figure 13 (Data at t - 51.5 min, RPLD model): Pinal fits are 

quite good for most of the experimental data. 

(iv) Figure 14 (Data at t ■ 51.5 min, IMGD model): Fits to the data 

are reasonably good. 

It should be kept in mind that n(r) data in Set 1 shows a kink in the 
curve around 2.0 pm, thereby indicating that the si«e distribution is 
bimodal. But since the analytic models used to fit the data arc unimodal, 

RPLD and IMGD models will not fit the data very well in the whole size range. 
Fitting with a two term bimodal analytic model was not attempted for Set 1 
n(r) data, even though it has been used in the retrieval of size distributions 
from T vs X data as discussed in Section 4. 3. 

It should further be noted that in Figs. 6-8, the data for r < 0. 3 pm 
are extrapolated data. No analytic fits were attempted for the Set 3 size 
distribution data for t ■ 2.5 and 62.0 min, shown in Fig. 9. 
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SIZE DISTRIBUTION 
SET 1 DATA 
TIME - 0.5 MIN 

MODEL 4 - IMGD 


O MEASURED DATA 

Q INITIAL ESTIMATES 
« 1.086E + 05 

P2 ■ 5.00 

P3 ■ 0.100 

P4 • 2.00 
O NLLS BEST FIT 
Iterations ■ 15 
p^ - 1.31E+02 
Pg - 5.75 

Pj ■ 0.100 

P4 ■ 1*93 


RADIUS r. pm 


FIGURE 12. Analytic representation of Set 1, t • 0.5 min size 
distribution data using an Inverse Modified Gamma 
Distribution (IMGD) Model. 
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SIZE DISTRIBUTION 
SET 1 DATA 
TIME - 51.25 MIN 

MODEL 2 - RPL 


O MEASURED DATA 

□ INITIAL ESTIMATES 
* 4.864E^06 

P2 - .290 

P3 * 8.00 

P 4 • 1.63 
O NLLS BEST FIT 
Iterations ■ 7 
p, - 5.633E + 06 
Pg ■ . 354 

P 3 - 5,72 
P 4 - 1.79 


-2 -I '> I 2 

RADIUS r, 

FIGURE 13. Analytic representation of Set 1, t - 51.25 min, size 

distribution data using an Regularized Power Law (RPL) Model. 
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SIZE DISTRIBUTION 
SET 1 DATA 
TIME - 51.25 MIN 

MODEL 4 - IMGD 


O MEASURED DATA 

□ INITIAL ESTIMATES 
- l.AeiE+.TS 

P 2 • 4.00 
Pj - 3.600E-02 
p^ « 3.00 
O NLLS BEST FIT 
Iterations ■ 7 

p^ ■ 1.08 

P 2 - 5.40 
Pg » 0.144 
P 4 - 2.16 


FIGURE 14. Analytic representation of Set i, t - 51.25 min, size 
distribution data us/ an Inverse Modified Gamma 
Distribution iXNGD) Model. 



4. RETRIEVAL OF AEROSOL SIZE DISTRIBUTION FROM 
MULT I SPECTRAL OPTICAL DEPTH MEASUREMENTS 


In Ref. 23 a method was described for retrieving sire distributions of 

aerosols from single wavelength maaeuremente of the decay of optical depth 

T (t) with time t. This method is valid under three assumptionsi (1) the 

particles are spherical and large (with radii r i V(P ** P pm, Ref. 24) i 

p m 

(2) they se .le in a quiet medium with a terminal velocity v^ given by Sto)ies' 
law (Ref. 25) t 


V 

s 



- Pp> 9 


( 1 ) 


where p^ and p^ are the specific gravities of the medium and particles, 

respectively, g is the acceleration due to gravity, and n is the medium 

m 

viscosity; and (3) scattering conditions are such that one can set Q “ 2.0, 
which is a reasonable assumption for large particles for which x > 6, where 
X “ 2‘lTr/A and Q is the efficiency factor. 

Any departures from the above conditions will yield errors in the results. 
For instance, if the particles are in stirred motion or undergoing microphysical 
<^anges (e.g. , coagulation, evaporation, growth, etc.), the extinction- 
sedimentation inversion technique of Ref. 23 is not applicable. In that case 
one must medce n independent measurements at any time t to retrieve n parameters 
needed to describe the size distribution. One simple way to achieve this is 
to ma)ce multispectral measurements of optical depth using n wavelengths. 

There are a number of ways of retrieving the aerosol size distribution 
from multispectral extinction (or optical depth) measurements, such as 
numerical Inversion techniques (e.g., Twomey-Chahine-type methods. Ref. 26), 

NLLS method of Box and Deepa)c (Ref. 27), and model dependent, table-search 
methods. The first two types of methods are iterative and can become quite 
time consuming; so, in the next three sections we shall describe fast 
approximate search methods. Comparisons of results will be made, wherever 
possible, with those obtained by the nonlinear least squares (NLLS) method. 

Teda.le Search Methods ; The optical depth T is defined by a Fredholm-type 
integral equation, namely. 



•2 *1 

where n(r) {cm Mm ] is the path-int^.grated sise distribution, Q is the 
Mie efficiency factor, r is radius (Mm ) , and m •* m* - im" is the complex 
aerosol refractive index. The method described here is similar to, but 
more general than, that of the search method of Box and Lo (Ref. 28) in 
which they used a modified gamma distribution (MGD) model (Ref. 21 ) , such 
as, Deirmendjiem's Haze I! model. In their method, they use a two-step 
approach — first step is to plot the T vs. X data on a semi-.iog or log-log 
graph paper and approximate it with a straight line, and then use a one- 
parameter, table-search method. But if the range of wavelengths is large, 
such a straight lino approximation is not always possible. For such a case, 
we developed the two-parameter and three -parameter taole-search methods 
(Refs. 29-30) which can deal with optical depth data for any range of wavelengths. 

In these methods, sets of tables for the T vs. X values were generated 
for the two or three adjustable parameters of the size distribution model 
and the best fit parameters for a particular experimental data set were 
found by searching the tabulated data base to determine the combination 
of parameters that yield the best fit. 

In addition, to the 2-parameter (Section 4.1) and 3-par2uneter (Section 4.2) 
methods, we extended the latter method to the case of two- term bimodal models, 
which is described in Section 4.3. 

In Section 4.4, we discuss the effects of noise on percentage errors in 
the size distributions retrieved by the direct Inversion (Ref. 23) of the 
single wavelength simulated data for T versus t for the case of aerosols 
settling quietly under gravity. 



4.1 FAST TABLE SEARCH (FTS) METHOD FOR RETRIEVAL OF AEROSOL 
SIZE DISTRIBUTION FROM MULTISPECTRAL EXTINCTION 
MEASUREMENTS! TWO- PARAMETER MODEL 


In the 2-parameter method, we wish to obtain from the expet imontal 
t vs. A data, two parameters that characterise the size distribution, such 
as, the total number of particles in the column a and the mode radius 
parameter b. To determine a and b we assume w analytic model for the aerosol 
size distribution H(r), with a and b as adjustable parameters, and then 
determine the combination of a and b by searching a precomputed tabulated 
data base of T(A) — generated for a judicially selected mesh of a and b values — 
to obtain the best fit to the experimental T(X) data. 

The tables should be generated for the same set of wavelengths with which 
the observations are made. Sensitivity of the search method retrievals to 
discrepancies in two iiets of wavelengths — the one to generate the tables and 
the other for measurements --is discussed by Box, Box and Deepak (Ref. 30) . In 
practice it was found easier to use the normalized path-integrated size distri- 
bution n (r) to calculate the tables, whore n (r) is defined here as that size 

° -2 ° 
distribution for which a ■> 1 pm , so that 

GO 

n (r) dr » 1 (3) 

o ° 

The values of i^(A), the normalized optical depth, are then calculated for a grid 
of values of b parameter of the selected n(r) model, using values of m * m' - im* 
corresponding to the wavelengths in the set. This procedure gives us a table 
of T^(A) for preselected values of b. 

The parameter a is defined by 

-2 

n(r) dr « a, [pm J (4) 



so that 


n(r) - an„(r) 

O 

— 3 —3 —1 

Note that n(r) (pm ] is related to the size distribution n(r) [cm pm ] as 


(5) 


n(r) ■ n(r) 10 /L(cm) 


jlO®/L(cm)J 


( 6 ) 


where L(cm) is the length of the aerosol chamber traversed by laser beams. 
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For a given value of b, using the norntalized optical depth tables one can 
calculate a from the relation 


a 


^i-1 


T<Ai> T^(X^) 


N 


^i-1 "^o ^^i^ 


(7) 


where T(A^) are the measurements taken for wavelength set A^ (i >■ 1 , 2 , . . . N) , 
and N is the number of wavelengths in the set. 

The best value of a is calculated for each value of b using Eq. (7). To 
determine which pair of a and b values give the beat fit to the T(A) data, 
we calculate X2» squares of the residuals, defined as 


N r 1' 

“I T(A.) - at’ (A.) 

i»l L ^ _ 


( 8 ) 


2 

for each pair of parameters, and determine which pair give the minimum x 

2 

value. The usual pattern is si that as b increases, X will decrease to some 

value, and then begin to increase. The best-fit a,b pair corresponds to that 
2 

point where X is a minimum. 


Examples of the Search Method 

The table-search method is illustrated by the use of two size distribution 

models, the Modified Gamma Distribution (MGD) and the Inverse Modified Gamma 

Distribution (IMGD) models (Ref. 21) . The MGD has an exponential fall off for 

large particles with radii r >> r^, the mode radius; whereas the IMGD has a power 

law fall off for r >> r . A brief description of the two models is given as 

ro 

follows. 


Modified Gamma Distribution (MGD) Model : 
In Ref. 21, the MGD model is defined as 


n(r) 


Pi' 


exp 


P4 

(-P3r ) 


( 9 ) 


where the parameters p and P . “ c are fixed; p is the mode radius parameter b 

-p 

(units, pm ^) ; and p^^, the scaling parameter, is chosen so that the total number 
of particles is a. For Eq. (9), the mode radius r^^^ is defined by 
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( 10 ) 




VP, 


and the kth moment* by 

"k ■ <Pl''P2’'’3 


-(P„ + k/P^) 


npj* + VP4) 


where 

pj4 ■ (Pjt i)/p^ 

From Eq. (10 ) , one obtains 


/ \ 
P 


"P. 


m 


and from Eq.(ll), the total number of particles a or p as 


a ■ M ■ 
o 


“P 


24 


r(P24> 


( 11 ) 


( 12 ) 


(13) 


(14) 


Then substituting for Pj^, p^ and p^ ■ c into Eq. 9 one can rewrite n(r) in the 
form 


n(r) ■ ac 


^P4 


r exp 




(15) 


where - (p^ + D/c. 


Then the effective radius is given by 


r “ b 
eff 


c 

l/C 

r 

P 2 + 4’ 

/ 

/ \ 
P 2 - 3 

N 

c 

. < 


c 

, 4 


(16) 


The special cases discussed here are obtained by giving different values to p^ 
p^, shown here as follows. 


a. Haze H: Pj “ 2, * 1 

, k 1 . 3 2 -br 

n(r) ■ b r e 

b. Haze L: p^ - 2, p^ - 0.5 

* » a .6 2 -b*^ 

n(r) = 240 ® 


(17) 


(18) 


and 
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( 19 ) 


c. Haze CHj - 2, - 3 


2 -br' 

n(r) ■ 3 a b r e. 


Inverse Modified Gamma Distribution (IMGD) Model s 
In Ref. 29, the IMGD model is defined aa 


r ^4\ ^2 

n(r) ■ exp|- p^/r J/r 


( 20 ) 


in which p^ and P2 “ c are fixed; p^ is the mode radius parameter b (units, Dm ^) 

and p. is chosen so that the total number of particles is a. The mode radius r 
•L m 

is given by 


1/P. 


*^m - <P3P4/P2> 


( 21 ) 


and the kth moment M^by 


-<P4j - k/p,) 


np,j - k/p^) 


( 22 ) 


where 


P 42 - 'Pj - 


(23) 


and a, the total number of particles, is given by 

■P. 


M 

M “ — 

o l_P4J 


42 


r(P42> 


Substituting for p^ in Eq. (20), Pj^ in Eq. (23), and p, 
rewrite n(r) as 


• P4 r 

/42 
c 1 

c 


'■^42 


P 4 



(24) 

c in Eq. (19), one can 

(25) 


The effective radius is then given by 


'e£f " *’1^) 


tiMhr] 


(26) 
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and the total mass concentration by 


M 

c 


p . b’ 
31, " 


[vp 


( 27 ) 


Prom Eqs. (25) and (26), one can easily see that if c £ 4, then Bnd 

cannot be defined. 


One special case used in the two-parameter method is for 
p^ ■■ 4 and p^ •> 1, which gives 


n(r) - ^(e"*’/’^)/r'* (28) 

-4 

v/ith r fall-off for r >> r . This will be referred to as H1NV41 model. 

m 

Computation 

In the formulation in Eqs. (7) and (8), it has been assumed that all 
measurements have been given equal weighting so far as the errors are con- 
cerned. (The case in which measurements are weighted according to the experi- 
mental err»jr will be considered in Section 4.2.) 

The tables of T^(^j^) were computed for five wavelengths given in Table 1. in Sec. 
and for the parameter values for MGD and IMGD models given in Table 3. 

The parameters a and b for each of four models (Haze H, Haze L, Haze CH and 
HINV41) were determined by the aforementioned method, and the results are given 
in Table 4. To determine how good these results are, numerical inversions of 
the data were also performed using a nonlinear least squares (NLLS) method and 
the results are also shown in Table 4. The approximate estimates for a and b 
were used as initial guesses for a and b in the NLLS program. 

The mode radii and the modal values of the size distributions at the 
different times are given in Tadiles 5 and 6, respectively. The size distributions 
obtained by us.lng the approximate estimates of the parameters for the different 
times during the experiment are plotted in Figs. 15 to 18, along with the experimental 
data for size distribution measurements ta)(en at times t ■ 0.5 and 51.25 min 
by means of a ten-stage Andersen Cascade Impactor. 
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TABLE 3. VALUES OP THE PARAMETERS p. FOR GENERATING THE TABLES FOR T(X) 


Parameters 

MGD 

INGD 

P2 

1.0, 2.0 and 3.0 


P4 


1.0, 2.0 and 3.0 

c 

0.3, 0.5, 1.0, 1.5 

3.5, 4.0, 5.0, 6.0, 7.0 


2.0 and 3.0 

and 8.0 

b 

0.1 (0.1) 1.0 (0.5) 5.0 

0.1 (0.1) 1.0 (0.5) 5.0 
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TABLE 4 . the RETRIEVED PARAMETERS a AND b BY SEARCH AND NLLS METHODS 


Time 


Search Method 



NLLS 


(min) 

a 

b 


a 

b 


Haze H 

0 

14.100 

15.00 

0.0154 

13.000 

14.60 

0.0017 

20 

3.130 

9.00 

0.0568 

3.020 

8.84 

0.0182 

40 

1.510 

7.00 

0.1840 

1.630 

7.23 

0.1102 

52 

1.400 

7.00 

0.1830 

1.340 

6.83 

0.1391 

80 

Haze L 

1.190 

7.00 

0.2750 

1.010 

6.41 

0.2256 

0 

40.100 

19.00 

0.0277 

39.600 

19.00 

0.0017 

20 

6.530 

14.00 

0.2870 

5.760 

13.60 

0.1927 

40 

3.820 

13.00 

0.5450 

3.030 

12.30 

0.4356 

52 

2.500 

12.00 

0.5340 

2.460 

11.90 

0.4624 

80 

2.120 

12.00 

0.5940 

1.790 

11.50 

0.5314 

Haze CH 

0 

5.700 

25.00 

0.0092 

5.760 

25.30 

0.8094 

20 

2.210 

10.00 

0.0086 

2.120 

9.40 

0.0071 

40 

1.280 

6.00 

0.0021 

1.280 

6.02 

0.0048 

52 

1.050 

5.00 

0.0015 

1.040 

4.92 

0.0066 

80 

0.773 

4.00 

0.0085 

0.704 

3.44 

0.0032 

HINV41 

0 

2110.000 

0.05 

1.5600 

2894.000 

0.0447 

0.9520 

20 

21.400 

0.25 

1.5600 

17.200 

0.2760 

1.2770 

40 

5.650 

0.40 

1.5700 

4.690 

0.4410 

1.3690 

52 

3.990 

0.45 

1.4300 

3.270 

0.5010 

1.2540 

80 

2.170 

0.55 

1.2600 

1.950 

0.5900 

1.1240 
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TABLE 5. MODE RADZZ OF RETRZEVED SZZE DZSTRZ BUT ZONE. 

The Experimental Values of Mode Radius at t ■ 0.5 and 51.25 min 
are 0.22 and 0.3 pm, respectively 


Model 

Inver- 

sion 

method 


Mode 

Radius (pm) 




t * 0 min 

20 min 

40 min 

52 min 

80 min 

Hase H 


Search 

0.1333 

0.2222 

0.2857 

0.2857 

0.2857 


NLLS 

0.1370 

0.2260 

0.2770 

0.29 U 

0.3120 

Haze L 


Search 

0.0443 

0.0816 

0.0947 

0.1111 

0.1111 


NLLS 

0.0445 

0.0865 

0.1060 

0.1120 

0.1210 

Haze CH 


Search 

0.2991 

0.4058 

0.4811 

0.5112 

0.5506 


NLLS 

0.3010 

0.4180 

0.4840 

0.5170 

0.5820 

HINV41 


Search 

0.0125 

0.0625 

0.1L90 

0.1125 

0.1375 


NLLS 

0.0112 

0.0691 

0.1100 

0.1250 

0.1480 



TABLE 6. NODAL VALUES OF SIZE DISTRIBUTIONS. 

Th« Experlmontal Value* of **<*'g,) ^ •«<* 51.25 ndn 

are 2.66(7)* and 5.0(6)* MR)/ raapactivaly 


Model 

Inver- 

sion 

method 



"'V 




t ■ 0 min 

20 min 

40 min 

52 min 

80 min 

Haze H 


Search 

4.887(7) 

6.516(6) 

2.451(6) 

2.269(6) 

1.935(6) 


NLLS 

4.380(7) 

6.180(6) 

2.720(6) 

2.120(6) 

1.500(6) 

Haze L 


Search 

2.415(8) 

2.138(7) 

1.078(7) 

6.004(6) 

5.108(6) 


MLLS 

2.371(8) 

1.778(7) 

7.680(6) 

5.856(6) 

3.936(6) 

Haze CH 


Search 

1.650(7) 

4.718(6) 

2.305('>) 

1.778(6) 

1.221(6) 


NLLS 

1.670(7) 

4.430(6) 

2.230(6) 

1.760(6) 

1.060(6) 

HINV41 


Search 

8.431(10) 

1.711(8) 

2.826(7) 

1.774(7) 

7.889(6) 


NLLS 

1.300(11) 

1.250(8) 

2.130(7) 

1.310(7) 

6.600(6) 


*The number within the parentheses refers to the power of 10, 
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ON(R)/OR 



TIME (MIN) 

O 0.0 
□ 20.0 
O 40.0 

A 52.0 
80.0 

TIMES FOR MEASURED 
SIZE DISTRIBUTION 
9 0.5 MIN 
S 51.25 MIN 


RADIUS, ym 


FIGURE IS. Retrieved size distributions for different times for 
Set 1 data using 2-parameter Haze H Model. 
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FIGURE 16. Retrieved size distributions for different times for 
Set 1 data using 2-parameter Haze L Model . 
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FIGURE 17, Retrieved size distributions for different times for 
Set 1 data using 2-parameter Haze CH Model. 
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Discussion of Rcsul ts 


I* Haze H Modal 

Hie retrieved size distribution results for the Haze H model are 

plotted in rig. IS, for the different times at which the experimental T(Xj^) data 

was available. This shows that as the time increases, r increases, whereas the 

n 

total number of particles decreases. Comparisons of the retrieved n(r) curve 
for 0.0 min with experimental data taken at t •> 0.5 min shows a fairly good 
agreement, the "true" mode radius is somewhat larger than retrieved value, but 
the modal values of n(r) are of the same magnitude. In the range 0.2 < r < 1.0 pm 
there seems to be a good agreement between the retrieved and measured results. 

Comparison of the retrieved curve at 52.0 min and the data taken at 
51.25 min is not so good. Their mode radii compare rather well, but their modal 
values differ, the true value being somewhat higher than the measured value. 

The rate of fall-off agrees reasonably well in the range 0.7 < r < 2.0 pm. 

In Table 5 , comparison of the search and NLLS results show quite a 
good agreement: b values differed by less than 5% in all cases, except for 

t ■ 80 nin, for which they differed by leas than 10%; 'a' values differed- by less 
than 10%, except for t » 80 min for which the difference was less than 20%. The 
higher discrepancies in a values are due tothe fact that in the search method 
the value of a depends on b, and hence errors in b will propagate and enhance the 
errors in a. Similarly, the mode radii differ by less than 5%, except for 
t •> 80 min for which they differed by 10%; and the modal values by less than 15%, 
except for t « 80 min, for which the difference is 30%. Comparisons with the 
experimental results show that both search and NLLS methods underestimate 
the mode radii and overestimate the modal values at t « 0.0 and 52 min. At 
52 min the mode radius estiiuates are good, but the maximum is underestimated by 
a factor of 2. 

2. Haze L Model 

The retrieved results for the Haze L model and experimental data points 
are plotted in Fig. 16. The shift in mode radius to larger values with increasing 
time and the corresponding drop in the total nuihber of particles is evident, but 
neither the 0.0 min nor the 52.0 min curves agree very well with the measured 
values. 

As was the case for Haze H, the agreement between the search and NLLS 
methods is quite good. The discrepancy between the b values was less than 6% 
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and the discrepancy between the a values was generally less than 20%. The 
discrepancies in the mode radii were less than 10%, but ranged up to 40% for 
the modal values. The calculated mode radii were much less than the "measured" 
mode radii for both methods for which data were available, but the modal values 
agreed quite well. 

3. Haze CH 

The results for this model are plotted in Fig. 17. The increase in mode 
radius with increasing time and the decrease In the total number of particles is 
evident. The agreement between the 0.0 min curve and the measured size distribu- 
tion at 0.5 min is quite good for 0.07 < r < 0.6 pm, although the mode radius is 
about 35% higher than the"true"mode radius and the maximum Is nearly 40% below 
the "true" maximum. The agreement between the 52 min and the measured values at 
51.25 min is not quite as good. The mode radius is 70% higher than the true 
mode radius and the maximum Is 20% higher than the true maximum. 

Comparison of the parameters obtained by the search method and the NLLS 
inversion shows that agreement is quite good, within 5% in most cases. Agreement 
between mode radii and maxima for the two methods is alsc^ generally, within aboul^ 
5%. 


4. H1NV41 

The results for this fit are plotted in Fig. 18. Tbis model shows the 
shift to larger mode radii with increasing time and also the drop in the number 
of particles, but agreement between the measured and calculated size distribu- 
tions is not good. The mode radii for the calculated distributions are higher 
than the true values and the maxima are much higher than the true maxima. 

Comparison of the results obtained by the search method and the NLLS 
inversion show that the estimates for b differ by about 10% while the estimates 
for bi' differ by up to 30%. The estimates of mode radii differ by about 10% 
while the differences in the maxima are of the order of 35%. These differences 
are much larger than w'^re encountered with the Modified Gamma Distributions. 

Conclusions 

The results presented here indicate that the two-parameter search method 
is a practical method for obtaining estimates of the mode radius parameter b and 
number of particles a for the simple models considered so far. The agreement with 
the results for numerical NLLS inversion method was quite good an'^. could easily 
be improved by using a smaller step-size for b in generating the tzdsles. 
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Another distinct advantage of the search method is that it is much faster 
than the NLLS method. For example, vting the search method, 96 sec of CP time 
was required to confute the tables and rhen determine the parameters for each 
of the five data sets used for four models. Using NLLS inversion, the times 
to determine the parameters for each of the five data sets for only one model 
ranged from 120 to 170 sec of CP time. 

However, in order to improve the accuraby of the retrievals by the search 
method, we decided to extend this technique to determine three adjustable 
parameters, so that the method can be extended to determine the fali-off rate 
for r » r^. The three-parameter method is described in the next section. 
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4.2 FAST TABLE SEARCH (FTS) ^<ETHOD FOR RETRIEVAL OF AEROSOL 
SIZE DISTRIBUTION FROM MULTI-SPECTRAL EXTINCTION 
MEASUREMENTS: THREE-PARAMETER MODEL 

In thii. section, we describe a three-parameter search method for deter- 
mining aerosol size distributions from multispectral optical depth T (X) 
measurements, by retrieving three parameters characterizing the size distri- 
bution. It is an extension of the two-parameter method described in 
Section 4.1. Here, we attempt to obtain three par^uneterB: the total number a 
of particles in the column, the mode radius parameter b, and the slope pareuneter 
c (corresponding to the rate of fall-off for large r) . To determine a, b, 

and c, we assume an analytic model for the path-integrated size distribution 
-2 *1 

n(r) [cm ym ], with a, b, and c as adjustable parameters, and then deter- 
mine the best combination of a, b, and c by searching precomputed sets of 
tables of T (X) for various: combinations of mesh of values of the three 
parameters to obtain the best fit to the experimental T(X) data. 

To calculate the tables of X (X) , we again work with the normalized 
size distribution n^(r), defined by Eq. (3). Similarly, the path-integrated 
size distribution n(r) is defined by Eq. (5). Thus, using n^j{r), tables for 
the normalized optical depth x^(X) cm be calculated for different values of 
b and c as given in Table 3. Each table is for one value of c and consists 

of X (X) values for different b values, 
o 

It will be assumed here that the measured optical thickness X(X^) is 
given by 


x(X^) - a T^(X^) + 


(29) 


where is the error. 

For a given value of b, the parameter a can be found from the relation 


^ 1-1 

rN - 2 , » . 2 

^i»l O 


where N is the total number of independent measurements and are the 
experimental errors in the measured values of X(X^). 


(30) 
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The best-fit parameter eet of a, b, and c for a given set of T (X) 

measurements is the combination which minimizes the weighted mean sum of 

2 

squares of the residuals X defined as 


111 w" 

M rx(x^) - 

ii.i L o, J 


( 31 ) 


2 

If X 5 the calculated values are consistent with measured values of the 
optical thickness, within the experimental error. 

After the best fit parameters have been obtained, one can easily calculate 

3 

the effective radius r (Vim) and total mass concentration N_(g/pm ) from 

ezf c 


eff 


M3/M2 


J^r^n(r)dr 

r^n(r)dr 


(32) 


M 

47Tp "3 
3 L 


(33) 


where 2 uid M. are the second euid third moments of the size distribution, 

2 3 

p io the particle density, and L is the beam path length within the aerosol 
media. 


Computations 

It was assumed that the error in the optical thickness measurements was 

2 

±3%. In calculating a and X « measurements were weighted according to 
the assumed experimental error. In Section 4.2 for the two-parameter, all 
measurements were given equal weighting resulting in a bias toward the small 
wavelengths, i.e., good estimates at X <» 3.39 pm, for which the T is much 
smaller. 

In generating the tables of T values, we assume analytic models for 
n(r) , such as. Modified Gamma Distribution (MGD) and Inverse Modified Gamma 
Distribution (IMGD) , and use pre-selected values for the parameters p^, b 
and c given in Table 3. The aerosol refractive values for four wavelengths 
(0.4416, 0.6328, 1.15, and 3.39 pm) are given in Table 2. 
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Measurements 


The measuranents of t(X, t) for Set 1. Set 2, and Set 3 data <Table 2) 
ara shown in Figs. 3-5, and the measurements of sixe distribution for Set 1 
and Set 3 are shown in Figs. 6 and 9. 

Discussion of Results 

The three-parameter search method was used to invert the T(X, t) data 
for Set 1, Set 2 and Set 3 in order to retrieve aerosol sixe distributions 
represented by two analytic models— the MGD and the INGD. The results for 
the three sets are discussed as follows. 

Set 1 Data - MGD Model : The results for the observed and calculated 

T, show that they are in reasonably good agreement. The results for the 

2 

goodness of fit for the MGD model are shown in Table 7, which gives X 

values for the three different values of the MGD parameter p^ at different 

times t «= 0, 20, 40, 52, and 80 min during the experiment. Remembering that 
2 

X ^ N indicates an acceptable fit and that N » 3 for the particular set of . 
data, we see that, in general, all the three values of p^, with p^ •* 2.0 
being slightly better than the other two. 

The measured and the calculated values of the mass concentrations, as 
functions of time are shown in Table 8. The calculated values are 15% to 
40% higher than measured values, with their difference increasing as time 
increases. Measurements seem to indicate that the aerosol sixe distribution 
is bimodal and this may contribute to over-estimation of the total mass 
concentration. The calculated values of M for the three values of p. differ 

C 4C 

by less than 2%. 

The effective radii of the time varying size distribution are shown in 

Table 9. The r values at any time for different p- values differed by 
eft 2 

less than 10%; but r Increases as time increases. The r values for 

eff eff 

the measured size distributions at t * 0.5 and 51.25 min differed by less 
than 5% from values calculated from retrieved results. 

The measured and retrieved size distributions are plotted in Figs. 19 
to 21 for three values of p,. All the three figures show that as time 
increases, the size distribution shifts toward larger particles and number 
density decreases. The exponential fall-off of MGD for larger r is too 
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TABLE 7. GOODNESS OF FIT FOR THREE VALUES OF MGD PARAMETER 
AND FOR SET 1 DATA 


Tima 

(min) 




?2 - 1.0 

Pj " 2.0 

Pj " 3.0 

0 

0.36 

0.27 

9.41 

20 

1.45 

2.11 

2.64 

40 

2.73 

5.85 

0.13 

52 

6.10 

0.10 

3.27 

80 

13.00 

2.56 

5.06 


TABLE 8. COMPARISON OF MEASURED AND CALCULATED MASS CONCENTRATIONS M 

C 

FOR SET 1 


Measured Calculated 


Time 

(min) 

g/m"* 

Time 

(min) 


g/m^ 

■ 

P2 - 1.0 

P2 - 2.0 

P2 - 3.0 

0.50 

0.653 

0 

0.761 

0.748 

0.758 

20.50 

0.606 

20 

0.734 

0.726 

0.732 

40.50 

0.541 

40 

0.681 

0.686 

0.681 

51.25 

0.492 

52 

0.688 

0.677 

0.687 
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TABLE 9. EFFECTIVE RADIUS AS A FUNCTION OF TIME FOR MGD AND SET 1 


Time 

(min) 


Effective 

Radius 


From 





measured 





size 





distribution 


Search Method 




?2 - 1.0 

Pj “ 2.0 

P2 = 3.0 

0 

0.366 

0. 365 

0.380 

0.357 

20 


0.510 

.0.507 

0.536 

40 


0.583 

0.602 

0.536 

52 

0.624 

0.656 

0.634 

0.655 

80 


0.656 

0.697 

0.714 
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FIGURE J9. Retrieved size distribution for different times 

Set I data using 3-parameter MGD Model with “ i-0; 
ar\d measured size distribution. 
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FIGURE 20. Retrioved size distributions for different times for 
Set 1 data using 3~parameter MOD Model with * ^>0; 
and measured size distribution. 
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FIGURE 21. Retrieved size distributions for different times for 
Set 1 data using 3-parameter MOD Model with “ 3.0; 
and measured size distribution. 
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steep to match with the measured distribution, but the MOD for the two cases 
of p^ ■ 2.0 and 3.0, does model the small radius (r 0) region well, 
especially for t •• 0.0 min case. 

CompMriMon of KoMultt Obtained with the Two~ and Three- Parameter 
Search Method* t In the three-parameter method, the best value for c was 3.0. 
In the two-parameter method (Sec. 4.1), three models were considered, one of 
which, namely Haze CH, had parameters p^ " 2.0 and p^ (or c) ■ 3.0 and is, 
therefore, similar to the three-parameter case with P2 " 2.0, as shown in 
Table 10. One can easily see in Table 10 that even though the two-parameter 
fit used equal weighting for all data points, yet the retrieved results are 
consistent with those retrieved using the three-parameter method, in which 
measurements were weighted according to the assumed experimental error of 
+ 3%. 


Set 1 Data - IMGD Model: The results for the observed and calculated 

2 

optical thickness show that they arc in poor agreement. Table 11 gives X 

2 

for the different values of p^ at the different times. In all cases, X ^ h 

indicating a poor fit, although there is some improvement as p^ increases and 

the fall-off of size distribution as r 0 becomes sharper. 

Table 12 gives the measured and calculated mass concentrations as a 

function of time. The calculated mass concentrations are consistently 25% 

to 40% higher than the measured values. The value of p^ appears to liave little 

effect on the calculated concentrations since the values are in agreement to 

within 2% except at 0.0 min where they differ by up to 10%. 

Table 13 lists effective radius as a function of time for the different 

values of p.. The value of p. has little effect on the effective radius; 

4 4 

the difference between the calculated radii is less than 10% except at 0.0 
min where it is 25%. The radii for p^ * 2.0 and p^ ■ 3.0 differ from the 
values for the measured size distribution by less than 10%. 

Figures 22 to 24 are plots of the measured and calculated size distri- 
butions for the different values of p^. In each case, the best fit for the 
calculated distributions was with c ■> 8.0. For p^ ■ 1.0 and p^ * 2.0, this 
was in good agreement with the measured distribution in the region r 
For p^ ■ 3.0, it was a little •steep. 
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TABLE 10, COMPARISON OF TWO-PARAMETER MODEL WITH THREE-PARAMETER 
MODEL FOR SET 1 


Time 

(min) 

Two- 

-parameter 

model 

Three 

-parameter 

model 

a 

r * 

m 


a 

b 

2 

X 

0 

5.700 

0.299 

0.009 

5.590 

0.30 

0.272 

20 

2.210 

0.406 

0.009 

2.290 

0.40 

0.110 

52 

1.050 

0.511 

0.002 

1.090 

0.50 

0.099 

82 

0.773 

0.551 

0.009 

0.757 

0.55 

0,560 


* r » mode radius. For the two-parameter model, this is given by 
in 

(2/3b)^'^^ where b is the mode radius parameter. 
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TABLE 11. GOODNESS OP PIT POR XMGD AND SET 1 


Time 

(min) 




P4 - 1-0 

P4 - 2.0 

P4 - 3.0 

0 

53.1 

19.8 

16.80 

20 

57.0 

20.5 

13.00 

40 

60.3 

22.8 

8.89 

52 

63.1 

16.5 

10.90 

80 

65.6 

29.0 

11.90 


TABLE 12. 

COMPARISON 

OP MEASURED AND CALCULATED MASS CONCENTRATIONS 


POR SET 1 





Measured 




Calculated 


Time 

(min) 

Mass 

cone. 

(g/m^) 

Time 

(min) 


Mass Concentration 


P 4 - 1-0 

P 4 - 2.0 p 

^ * 3-0 

4 

0.50 

0.653 

0 

0.839 

0.763 

0.785 

20.50 

0.606 

20 

0.741 

0.728 

0.728 

40.50 

0.541 

40 

0.683 

0.672 

0.678 

51.25 

0.492 

52 

0.697 

0.685 

0.690 


52 


TABLE 

1 3. EFFECTIVE 

RADIUS AS A FUNCTION 

OF TIME FOR IMGD 

AND SET 1 



Effective 

radius r^^^(pm) 


Time 

From 

measured 


Search method 


(min) 

size 

distribution 

P 4 * 1-0 

P4 “ 2.0 

P4 - 3-0 

0 

0.366 

0.30 

0.376 

0.343 

20 


0.50 

0.527 

0.549 

40 


0.60 

0.602 

0.617 

52 

0.624 

0.70 

0.677 

0.686 

80 


0.70 

0.752 

0.686 
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Retrieved distributions for different times for 
Set 1 data using 3-parameter IMGD Model using 3 

and measured size distribution. 


The best ao:<recinont between the measured and calculated distribuitons 
at 0.5 min was obtained with p- ■ 2.0 and at 51.25 min with p. ■ 1.0. The 
calculated distributions for p^ ■ 3.’0 were much narrower than the measured 
distributions. 

Set 2 Data - MOD Model : The results of applying the search method to 

this data are not so good. Table 14 gives the x values for the different 
values of p^ at different times. We see from this table that no value of 
p^ can be considered to be any better than any other, although at each time 
at least one distribution satisfies the criteria for an acceptable fit, 

^ 3. 

Figures 25 to 27 are plots of the calculated size distribution for the 
different values of p^. There is a shift to larger particles and a decrease 
in the total number of particles as the time increases. 

Set 2 Data - IMGD Model i The results of applying the search method to 

2 

this data are generally poor as is evident from Table 15 which gives X 

as a function of p^ and time. There is some improvement as p^ increases and 

the distribution becomes narrower but even p^ = 3.0 is a very poor fit. 

Figures 28-30 are plots of the calculated size distribution. All of 
the plots show a shift toward larger particles and a decrease in the total 
number of particles as the experiment progresses. 

Set 3 Data - MOD Model : The results of applying the search method to 
these data are in reason^lbly good agreement. In many cases, the observed 
and calculated optical thicknesses are in agreement to within the experi- 
mental error. There were four measurements in this data set and we see 

2 

from Table 16 that on a number of. occasions X ^4 indicating a poor fit. 

Comparison of the measured and calculated mass concentrations shows 
that, in most cases, they are in agreement to within 5% and the maximum 
difference is 12% (Table 17) . The value of the fixed parameter P2 has little 
effect on the calculated concentrations since they are in agreement to 
within 5%. This indicates that the Modified Gamma Distribution is a suitable 
model for obtaining the mass concentration for this particular data set. 
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TABLE 14. GOODNESS OF FIT FOR MGD FOR SET 2 


Time 

(min) 




P 2 ■ 1*0 

P 2 - 2.0 

P 2 - 3.0 

0 

15.40 

6.360 

0.393 

20 

4.26 

0.258 

11.300 

30 

13.00 

3.580 

2.690 

40 

1.58 

4.560 

0.334 

50 

4.58 

0.365 

4.030 
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TABLE 15. G(X)DNESS OF FIT FOR IMGD AND SET 2 


Time 

(min) 




P 4 - 1-0 

* 2.0 

P 4 - 3-0 

0 

62.4 

18.6 

8.16 

20 

57.7 

21.2 

19.80 

30 

68.3 

28.6 

16.30 

40 

55.0 

19.3 

6.87 

50 

56.7 

13.0 

8.78 
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FIGURE 29-. Retrieved distributions for different times for 

Set 2 data using 3-parameter JMGD Model with ■ 2.0 
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TABLE 16. GOODNESS OF FIT FOR MGD AND SET 3 


Time 

(min) 


2 

X 


?2 - 1-0 

P 2 “ 2-0 

P 2 " ^*0 

2 

25.80 

16.50 

14.300 

30 

15.30 

3.17 

14.000 

62 

3.75 

7.29 

0.697 

100 

9.23 

4.77 

4.600 

150 

6.54 

3.68 

3.330 
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TABLE 17. COMPARISON OF MEASURED AND CALCULATED MASS CONCENTRATIONS 


FOR SET 3 


Measured 




Calculated 


Time 

(min) 

Mass 

cone. 

Time 

(min) 


Mass concentration 
(g/m^) 

2 “ 3.0 

?2 = 1.0 

“2.0 p 

C.5 

0.76 

2 

0.733 

0.723 

0.738 

30.5 

0.72 

30 

0.736 

0.730 

0.747 

58.5 

0.65 







62 

0.636 

0.622 

0.638 

67.5 

0.62 





87.5 

0.59 







100 

0.548 

0.566 

0.558 

122.5 

0.46 





152.5 

0.40 

150 

0.436 

0.445 

0.442 
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The effective radii as a function of time for the different values oi 
the fixed parameter are given in Table 18. The radii increase steadily 
with time and the values differ by less than 10% for the different values 
of The effective radii at 2.5 min and 62.0 could not be determined 

from the sma'. 1 number of size distribution measureme.r .s to compare with the 
search method results. 

The measured and calculated size distributions are plotted in Figs. 

31 to 33. As time passes, there is a shift toward larger particles and 
the total number of particles decreases. For all values of calcu- 

lated distributions drop off faster than the measured size distributions 
as r Lack of measurements for r < 3 ym means that no conclusions can 

be drawn about how well the calculated distributions model Uie overall 
features of the aerosol size distribution. 

get 3 Data - JNGD Model : The results of using the search method to 

2 ^ 

fit these data are very poor. Remembering that for this data set X * 4 
indicates a good fit, we see From Table 19 that the fits are generally very 
poor. The measured and calculated mass concentrations are in agreement to 
within 10%, except at 150 min where the difference is 18.5% (Table 20). 
Table 21 shows the calculated and measured results of r 

ef f 

The observed and calculated size distributions are plotted in Figs. 

34 to 36. The calculated distributions show a steady shift toward larger 
particles and a decrease in the total number of particles with time. This 
is consistent with the behavior of the measured distributions but the 
calculated distributions fall off more rapidly as r 

Comparison of the Search Method 
and Nonlinear Least Squares 

The Modified Gamma Distribution was fitted to the T versus X data 
using a Nonlinear Least Squares (NLLS) code. The initial estimates for 
the parameters in the NLLS program were the best fit estimates obtained 
by the search method. The purpose was to see how close the NLLS final 
estimates were to the search method estimates. 

The formula used in the NLLS code was that given in Eq. 9 of . 

Section 4.1 since this is the simplest expression for the Modified Gamma 
Distribution. The parameters a and b can be calculated from Pj^, p^, p^* 
and p^ using Eqs. 14 and 10 of Section 4.1. 



TABLE 18. EFFECTIVE RALIUS AS A FUNCTION OF TIME FOR MGD AND SET 3 


Time 

(min) 


Effective radius (pm) 


P2 - 1.0 

" 2.0 

P 2 ■= 3-0 

2 

0. 376 

0.380 

0.407 

30 

0.583 

0.571 

0.595 

62 

0.656 

0.634 

0.655 

100 

0.729 

0.761 

0.746 

150 

0.802 

0.824 

0.814 
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FIGURE 31. Retrieved distributions for different times for 

:^et 3 data using 3'‘parameter MGD Model with * 1.0; 
and measured size distribution. 


70 


tio/ciiiNa 




OWCR}/Oft 



FIGURE 33. Retrieved distributions for different times for 

Set 3 data using 3-parameter MGD Model with p_ » 3.0; 
and measured size distribution. 
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TABLE 19. GOODNESS OP FIT FU« IMGD AND SET 3 


Time 

(min) 




P4 “ 1-0 

P 4 - 2.0 

P 4 “ 3.0 

2 

139.0 

74.6 

31.100 

30 

89.2 

35.4 

31.000 

62 

58.2 

22.2 

15.000 

100 

47.0 

10.9 

1.550 

150 

37.1 

4.37 

0.613 


TABLE 20. 

COMPARISON OF MEASURED AND CALCULATED MASS CONCENTRATIONS 
FOR SET 3 

Measured 



Calculated 


Time 

(min) 

Mass 

cone. 

(g/m^) 

Time 

(min) 

Mass concentration 
(g/m^) 

p^ - 1.0 

P 4 - 2.0 

P 4 “ 3.0 

0.5 

0.76 

2 

0.728 

0.710 

0.724 

30.5 

0.72 

30 

0.725 

0.737 

1 699 

58.5 

0.65 







62 

0.649 

0.644 

0.656 

67.5 

0.62 





87.5 

0.59 







100 

0.579 

0.562 

0.568 

122.5 

0.46 





152.5 

0.40 

150 

0.473 

0.450 

■ 0.449 
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TABLE 21 . EFFECTIVE RADIUS AS A FUNCTION OP TIME FOR IMGD AND SET 3 


TilM 

(min) 


Effectiv* radius (pm) 


P4 - 1-0 

- 2.0 

P 4 - 3.0 

2 

0.40 

0.376 

0.412 

30 

0.60 

0.602 

0.549 

62 

0.70 

0.677 

0.686 

100 

O.BO 

0.752 

0.754 

150 

0.90 

0.827 

0.823 
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FIGURE 35. Retrieved distributions for different times for 

Set 3 data using 3-parameter JMGD Model with •= 2.0 ; 
and measured size distribution. 
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FIGURE 36. Retrieved distributions for different times for 

Set 3 data using 3-parameter IMCD Model with * 3.0 
and measured size distribution. 
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The results Cor Set 1 are given in Tables 22 and 23, for Set 2 in 

Ted}les 24 «nd 25, end for Set 3 in Tables 26 and 27. Tot all throe sets, 

2 

the NU,S results for • 3.0 with X 3 and ingirobable values for 
ci'^ 20) BO these results will not be considered in the following discussion. 
The poor results for •> 3.0 may have been because the NLLS code was not 
suitable for fitting such a narrow size distribution. 

The NLLS estimates of ”a” for the three days were generally within 
10% of the search method estimates and the largest difference was 15%. 

For "b" the NLLS estimates were generally within 5% of the search 
method estimates although differences up to 15% were noted. The 
differences between the NLLS estimates of ”b'' and the search method 
estimates were less than the step size used in the search method except 
at 2 min on Set 3 where the NLLS estimate for "b“ with p^ - 1.0 is 
much higher. 

The biggest differences between the NLLS estimates and the search method 
are for c, where the NLLS cstimatof' > ; s consistently 20% to 30% lower than 
the search method estimates. The exception to this is in the case of t => 2 min 
on Set 3 for p^ ■* 1.0 where the NLLS estimate is double the search method 
estimate. The differences between the NLLS estimates and the search method 
estimates were generally less than the step size in c. 

The effective radii calculated using the NLLS estimates were 
generally within 10% of the values calculated using the search method 
estimates. For Set 1 the NLLS estimates of effective radius were 
up to 10% higher than the values calculated from the measured size 
distribution, compared to 5% for the search method. 

The mass concentrations calculated from the NILS estimates differed 
by less than 10% from those calculated using the search method except 
for Set 1 where it rose to 35% for p^ “ 2.0. The NLLS estimates of 
mass concentration were up to 50% higher than the measured mass concentra- 
tions for Pj « 1.0. This dropped to 30% for “ 2.0. 

The NLLS estimates of mass concentration for Set 3 differed from 
che measured values by less than 10% compared with 12% for the search 
method . 

The search method is much faster and thus cheaper to use than the 
NLLS. For both MGD and IMGD, the search method took 16 seconds to do 
45 complete table searches and generate the required reports. The 
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TABLE 22. COMPARISON OF SEARCH METHOD AND NONLINEAR LEAST SQUARES 


FOR SET 1 


Time Search method NLLS 


(min) 




2 





"2 


a 

b 

c 

X 

a 

b 

c 

X 

P2 “ 1. 

0 








0 

7.370 

0.24 

3.0 

0.356 

6.740 

0.233 

2.32 

18.370 

20 

2.590 

0. 35 

3.0 

1.450 

2.464 

0. 331 

2.44 

0.916 

40 

1.610 

0.40 

3.0 

2.730 

1.516 

0.380 

2.43 

0.723 

52 

1.140 

0.45 

3.0 

6.100 

1.159 

0.441 

2.76 

0.301 

80 

1.000 

0.45 

3.0 

13.000 

0.904 

0.443 

2.53 

0.657 

P2 - 2. 

0 








0 

5.590 

0.30 

3.0 

0.272 

5.144 

0.297 

2.45 

24.220 

20 

2.290 

0.40 

3.0 

2.110 

2.035 

0.401 

2.41 

3.338 

40 

1.440 

0.40 

2.0 

5.850 

1.407 

0.398 

1.83 

1.098 

52 

1.090 

0.50 

3.0 

0.099 

1.013 

0.494 

2.42 

0.050 

80 

0.757 

0.55 

3.0 

2.560 

0.740 

0.546 

2.64 

1.928 

P2 - 3. 

0 








0 

6.380 

0.30 

3.0 

9.410 

0.157 

1.069 

26.30 

1084.000 

20 

1.830 

0.45 

3.0 

2.640 

0.199 

1.074 

18.90 

803.000 

40 

1.240 

0.50 

3.0 

0.126 

0.233 

1.050 

27.90 

299.000 

52 

0.938 

0.55 

3.0 

3.270 

0.209 

1.090 

18.30 

405.000 

80 

0.664 

0.60 

3.0 

5.060 

0.193 

1.096 

18.00 

298.000 


4 
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TABLE 23. COMPARISON OF SEARCH METHOD AND NONLINEAR LEAST SQUARES 
FOR SET 1 


Time 

(min) 

Search 

method 


NLLS 


*^eff 

M 

c 

*^eff 


M 

c 

?2 ■ 1-0 

0 

0.3€5 

0.761 

0.3 '4 


0.813 

20 

0.510 

0.734 

0.543 


0.791 

40 

0.583 

0.681 

0.624 


0.738 

52 

0.656 

0.688 

0.672 


0.733 

80 

0.656 

0.605 

0.710 


0.657 

P, - 2.0 

0 

0.380 

0.748 

0.407 


0.797 

20 

0.507 

0.726 

0.553 


0.785 

40 

0.602 

0.686 

0.628 


0.683 

52 

0.634 

0.677 

0.680 


0.589 

80 

0.697 

0.624 

0.725 


0.463 

Pj - 3.0 

0 

0.357 

0.758 

0.987 


0.450 

20 

0.536 

0.732 

1.008 


0.606 

40 

0.595 

0.681 

0.966 


0.628 

52 

0.655 

0.687 

1.025 


0.668 

80 

0.714 

0.632 

1.032 


0.628 
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TABLE 24. CX>MPAR1S0N OP SEARCH METHOD AND NONLINEAR LEAST SQUARES 


FOR SET 2 


Time 

(min) 

Search Method 



NLLS 



a 

b 

c 

2 

X 

a 

b 

c 


Pp " 1 

.0 








0 

6.750 

0.25 

3.0 

15.40 

7.458 

0.241 

2.79 

3.508 

20 

2.530 

0.35 

3.0 

4.26 

2.532 

0.334 

2.58 

3.358 

30 

2.130 

0.30 

2.0 

13.00 

2.187 

0.295 

1.91 

0.496 

40 

1.640 

0.40 

3.0 

1.58 

1.541 

0.379 

2.42 

0.196 

50 

1.170 

0.45 

3.0 

4.58 

1.178 

0.438 

2.70 

0.032 

P2 - 2, 

.0 








0 

9.590 

0.25 

3.0 

6.36 

7.034 

0.250 

2.14 

33.376 

20 

2.230 

0.40 

3.0 

0.258 

2.111 

0.395 

2.48 

6.269 

30 

1.570 

0.45 

3.0 

3.58 

1.558 

0.445 

2.63 

3.929 

40 

1.460 

0.40 

2.0 

4.56 

1.430 

0.398 

1.83 

0.427 

50 

1.120 

0.50 

3.0 

0.365 

1.068 

0.499 

2.39 

0.396 

P2 - 3. 

0 








0 

7,730 

0.25 

2.0 

0.393 

0.178 

0.994 

21.50 

986.000 

20 

1.770 

0.45 

3.0 

11.300 

0.179 

1.088 

21.40 

853.000 

30 

1.830 

0.40 

2.0 

2.690 

0.197 

1.084 

17.70 

699.000 

40 

1.260 

0.50 

3.0 

0.334 

0.206 

1.087 

18.50 

543.000 

50 

0.959 

0.55 

3.0 

4.030 

0.216 

1.090 

18.30 

394.000 
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TABLE 26. COMPARISON OF SEARCH METHOD AND NONLINEAR LEAST SQUARES 
FOR SET 2 


Time 

(min) 

Search 

method 

NLLS 


‘^eff 

M 

c 

^eff 

M 

c 

Pj - 1-0 





0 

0.365 

0.697 

0.365 

0.758 

20 

0.510 

0.715 

0.530 

0.768 

30 

0.564 

0.714 

0.575 

0.761 

40 

0.583 

0.691 

0.625 

0.752 

50 

0.656 

0.705 

0.675 

0.751 

- 2.0 





0 

0.317 

0.742 

0.363 

0.754 

20 

0.507 

0.708 

0.538 

0.765 

30 

0.571 

0.710 

0.592 

0.764 

40 

0.602 

0.696 

0.627 

0.749 

50 

0.634 

0.693 

0.691 

0.812 

Pj “ 3.0 





0 

0.339 

0.721 

0.926 

0.420 

20 

0.536 

0.711 

1.014 

0.556 

30 

0.543 

0.698 

1.022 

0.624 

40 

0.595 

0.691 

1.022 

0.651 

50 

0.655 

0.702 

1.025 

0.689 
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TABLE 26. COMPARISON OF SEARCH METHOD AND NONLINEAR LEAST SQUARES 


FOR SET 3 


Time 

(min) 

Search method 



NLLS 



a 

b 

c 


a 

b 

c 


- 1 

.0 








2 

7.380 

0.20 

2.0 

25.800 

4.820 

0.325 

3.98 

7.853 

30 

1.740 

0.40 

3.0 

15.300 

1.769 

0.394 

2.79 

4.927 

62 

1.060 

0.45 

3.0 

3.750 

1.007 

0.431 

2.51 

4.213 

100 

0.664 

0.50 

3.0 

9.230 

0.625 

0.477 

2.47 

13.454 

150 

0.397 

0.55 

3.0 

6.540 

0.404 

0.484 

2.25 

13.473 

P2 - 2 

.0 








2 

5.400 

0.30 

3.0 

16.500 

4.037 

0.300 

2.35 

5.017 

62 

1.010 

0.50 

3.0 

7.290 

0.884 

0.494 

2.27 

1.160 

100 

0.529 

0.60 

3.0 

4.770 

0.520 

0.597 

2.73 

4.293 

150 

0.327 

0.65 

3.0 

3.680 

0.326 

0.648 

2.80 

3.575 

P2 » 3 

.0 








2 

4.580 

0.30 

2.0 

14 . 300 

0.181 

1.060 

26.30 

1661.000 

30 

1.360 

0.50 

3.0 

14.000 

0.221 

1.088 

18.30 

904.000 

62 

0.871 

0.55 

3.0 

0.697 

0.208 

1.095 

19.10 

431.000 

100 

0.562 

0.55 

2.0 

4.600 

0.173 

1.112 

11.70 

383.000 

150 

0.342 

0.60 

2.0 

3.330 

0.133 

1.119 

10.40 

269.000 
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TABLE 27. COMPARISON OF SEARCH METHOD AND NONLINEAR LEAST SQUARES 
FOR SET 3 


Tiiiie 

(min) 

Search 

method 

NLLS 


*^eff 

M 

c 

'eff 

M 

c 

Pj - 1.0 





2 

0.400 

0 •’28 

0.418 

0.727 

30 

0.600 

0. /25 

0.597 

0.785 

62 

0.700 

0.649 

0.695 

0.683 

100 

0.800 

0.579 

0.776 

0.589 

150 

0.900 

0.473 

0.837 

0.463 

Pj - 2.0 





2 

0.376 

0.710 

0.418 

0.678 

30 

0.602 

0.737 

0.601 

0.782 

62 

0.677 

0.644 

0.700 

0.688 

100 

0.752 

0.562 

0.783 

0.594 

150 

0.827 

0.450 

0.842 

0.465 

?2 - 3.0 





2 

0.412 

0.724 

0.978 

0.505 

30 

0.549 

0.699 

1.023 

0.702 

62 

0.686 

0.656 

1.028 

0.670 

100 

0.754 

0.568 

1.079 

0.639 

150 

0.823 

0.449 

1.096 

0.513 
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45 searches were for five sets of data for each of the three sets for 
each of the three values of the fixed parameters. 

The NLLS code is an iterative or? and thus the amount of time taken 
to determine the best fit parameters depends on how quickly it 
converges. A total of 48 fits were done. These correspond to the 45 
searches done by the NLLS, with three being repeated. The total time 
required for these fits was 1705 seconds which averages to about 
' 36 seconds per fit. As has been pointed out above the length of time 
taken for one NLLS fit depends on the number of iterations. The total 
number of iterations in the above 48 fits was 387 which is approximately 
4.4 seconds per iteration. Since none of the NLLS fits converged in 
less than three iterations, even with good initial estimates, the NLLS 
is much slower than the search method. 

Discussion and Conclusions 

The results presented in this report Indicate that useful information 
about aerosol size distributions ceui be obtained from multispectral measure- 
ments of optical thickness using the three parameter search method. It was 
found that the Modified Gainma Distribution generally gave a better fit to 
the T vs A curve and better estimates of the effective radius and mass 
concentration than the Inverse Modified G;anma Distribution. On the 
other hand, the Inverse Modified Gamma Distribution gave better agree- 
ment with the measured and c Iculated size distributions in the region r 

It is not clear from the results presented here whether the Inverse 
Modified Gamma Distribution is am unsuitable model or whether the problem 
lies mainly in the values selected for the parameters, particularly 
the fixed parameter p^. More work will be necessary in order to answer 
this question. 

The similarity between .the NLLS estimates and the search method 
estimates indicates that the search method is capable of giving good 
estimates of the parameters if the grids are well chosen. The fact that 
the NLLS estimates for c were consistently lower than the search 
method estimates '.ndicates that an extra value of c is needed, around 
2,5, in the search method grid. 
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The search method has a definite cost advantage over NLLS and the 
fact that it is much faster, and is not an iterative process, would be 
important if information about the size distribution were to be obtained 
from the T versus X measurements on an on-line basis during the course of 
an experiment. The search method estimates also provide good initial estimates 
for the NlXS. This is important for rapid convergence. 

The Modified Gamma Distribution generally falls off faster than the 
measured distributions as r “ but it seems liKaly that this will not 
be of great importance if the distribution gives good estimates of the 
effective radius and total mass concentration. This applies particularly 
when only a few measurements are avail 2 ible, as was the case here, since 
good estimates of these features indicate that the essential nature of 
the distribution has been captured and the retrieved distribution can 
be considered an equivalent distribution. 

At present, in determining the best fit parameters for the 3-parameter 
search method, the complete table for each value of c is searched, which is 
probably not the most efficient way. Further worh needs to be done on this 
aspect in order to make the three parameter method more efficient. 
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4.3 FAST TABLE SEARCH (FTS) METHOD FOR RETRIEVAL OF 
AEROSOL SIZE DISTRIBUTION FROM MULTI -SPECTRAL 
EXTINCTICM^ MEASUREMENTS: 

TWO-TERI- BIMODAL MODEL 


Since the measured size distribution in Fig. 6 for Set 1 shows a kink 
around r « 1 pm, thereby suggesting a bimodal distribution, it v.'as decided 
to see whether any improvement in fits to the T(X,t) data could be obtained 
by using a two-term analytic model representing a bimodal size distribution. 
The model chosen was the sum of two Haze H distributions, i.e.. 


n(r) 




2 -*“l' 
r e 


2 2 


2 

r e 


(34) 


where each term represents a mode of a bimodal size distribution, 2 md adjustable 
parameters are as follows: 

“ total number of particles in the first mode 
b^ « mode radius for the first mode 
a^ “ total number of particles in the^ second mode 
” mode radius for the second mode. 

Since data for only three wavelengths 'were available for Set 1 and 
there are four parvneters to be determined, it was decided to fix b^^ and b^ 
eutd then determine a^ and a^ from the data. If the distribution changes with 
time, this would presumably be reflected in changes in the relative proportions 
of the two modes. 

The values chosen for the b's were 15.0 for b^ and 7.0 for b^, which 
corresponded to the best fit values for Haze H fitted to extinction data 
taken at times 0 and 52 min. The values of a^ ^uld a^ were found by solving 
the three simultaneous equations 

T(X^) = *1^1 + *2^2 i “ 1. 2, 3 (35) 


in the least squares sense; where T(X^) are measured data, Tj^(X^^) and 7 
are the extinction values for the first and second modes, calculated by using 
the normalized size distribution n (r) , defined as in Eq. 3 in Section 4.1. 
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Tables of the T (A.) were calculated as described in Sections 4.1 and 4.2. 
o X 


The 


expressions for a^ and are then given by 


•l " 


I T,=*(X^) - E T,(X^)T,(X,) 

p(X^) tj^(X^) ?: Tj^(X^) ~ f ~^2^^i^ I 

I E Tj(X.) Tj(X^) 1^ - E Tj^Xj) E t/(X^) 


Equation (34) was also solved for a^ and a^ using a nonlinear least 
squares code. Since the system is in fact linear if neither bj^ nor b 2 are to 
be determined, this is not the most efficient method of solving the equation. 


Results 

The distributions obtained for 0.0 and 52.0 for both methods are given 
in Figs. 37 and 38 and the parameter estimates are given in Table 27. Figure 37 
shows that both methods give distributions which agree quite well with the 
measured distribution, especially for r > 0.5 pm. The mode radii of the 
calculated distributions tend to be lower than the mode radius for the data 
while the maxima are higher. 

Table 28 gives a comparison of the calculated and true extinction values 
for the different wavelengths. It can be seen that both methods give good 
estimates tor the extinction at 0.6328 pm and 1.15 pm but the nonlinear least 
squares estimate is better for the 3.39 pm data. It is interesting to note 
the differences in the values of a^^ and a^ obtained by the two methods. 

Figure 38 shows that the linear method gives a reasonedsly good fit to the 
measured size distribution while the nonlinear least squares (NLLS) fit is 
very poor. One problem with the NLLS fit is that all the values are a factor 
of 10 to 15 too low and also the run had not reached convergence in 10 
iterations. 

A problem to be noted in the case of both these methods is that for the 
52 min data, a^ is negati\e and this gives rise to nonphysical negative size 
distributions for r < 0.2 pm. The problem of getting negative parameters is 
always a possibility wherever two analytic terms are added as was done here. 
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TABLE 27 s PARAMETERS FOR SET 1 DATA AND TWO TERM HAZE H HAZE H MODEL 


Time 

(min) 


Search 

Method 

• 

NLLS 

0.0 

*1 

13.72 


2.281 X 10^ 


“2 

0.066 


18.0 

52.0 


-1.199 


-121.0 


*2 

I.'-' 7 


15.9 


TABLE 

28: COMPARISON OF T (X) VALUES 

FOR SET 1 DATA 





T(X) 


Time 

(min) 

X (ym) 

Measured * 

Search 

Method 

NLLS 

0.0 

0.6328 

6.8 

6. 81 

6.80 


1.15 

3.7 

3.67 

3.70 


3.39 

0.45 

0.53 

0.44 

52.0 

0.6328 

6.8 

2.85 

0.17 


1.15 

3.7 

3.47 

0.20 


3.39 

0.45 

0.99 

0.05 

*at times 

t s 0.0 and 

51.25 min 
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n{r) 



FIGURE 37. Size distributions for Set 1 data: ® denotes measurements 

at t “ 0.5 min; dashed line, retrievals by NLLS code 
at t ■ 0; and solid line, retrievals by search method 
at t “ 0. 
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0 SET I DATA 

TIW • i1.?5 mm 

SO MOOEl : 

HAU H 4 ha;e h 

so PAlUVIITtRS 

It I 4 52.0 vjm 

SEARCH method 
I, . -I.19J 
b, ■ TS.O 
Ij • 1.687 
bj • 7.0 

KLIS HETMOO 
ITERATIONS • 10 

ij . -121.0 

b, > 1S.0 
•g • 15.9 
bg - 7.0 


FIGURE 38. Size distributions for Set 1 data: 9 denotes measurements 

at t * 51.25 min; dashed line, retrievals by NLLS code 
at t * 52 min; and solid line, retrievals by search meth.. ’ 
at t * 52 min. 
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2 3 

Figure* 39 and 40 give the plot* of dN/d log d, d dN/d log d and d dN/d log d 
for 0.0 and 52.0 min, xespectively . When these plots are con^red with the 
corresponding plots for the measure J n(r) data given in Pigs. 7 and B it is 
seen that even with the suni of two Haze H terms, the bimodality in retrieved 
results is not too clear. 

The calculated curves for the 0.0 min data are a little broader than the 
measured curves but do peak at about the same diameters and have very «<imilar 
maxima. The calculated curves are very similar to the measured curves in 
the region 0.5 < d < 2 ym, but the calculated curves drop off too steeply 
beyond this point. 

The calculated curves for the 52.0 min data peak a higher diameter than 
the measured curves and although the maximum for the dN/d log d curve is lower 
than the measured value, the maxima for the other two curves are higher than the 
calculated values. 

Conclusions 

The use of a two-term model does appear to have some merit in determining 
aerosol size distribution from T(X) data, but it should be clearly stated that 
three independent measurements are really not enough to obtain accurate results. 
Certainly if details such as the bimodal peaks and valleys in the curves are to 
be determined, more wavelengths are necessary. However, for the sake of complete- 
ness we decided to include this case in the final report. 

Thera are a number of problems inherent in the method used here. Since 
only three wavelengths were available, it was necessary to fix the mode 
radius parameters for the two terms. To be done most efectively, this 
requires some a priori knowledge of the aerosol size distribution. 

Another problem is the possibility that terms will be subtracted, rather 
than added, as happened with the data for 52 min. This opens up the possi- 
bility of negative size distributions which is clearly a physical impossibility. 
Further work using greater number of independent measurements taken for different 
wavelengths, needs to be done to draw firm conclusions about the accuracy of 
the results. 
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so NODCL : 
ha;i n • H/.:i h 

TINt t • 0.0 MIN 

StARCH niTHOO 
• , • 13 .)? 
b, • 1S.0 
• 0.066 
bj • 7.0 


DIAMCTER d. pm 


FIGURE 39. Three representations (or moments) tor Fig. 37 
rettieved size distribution: number-, area-, 

and volume-, log size distributions. 
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SB MOra: 

H42E N * HAIt H 

TIK t • 52.0 

SUUKH HTTHOO 

a, • *1.199 

b, > 15.0 
a^ • 1.517 
bj . 7.0 


FIGURE 40-. TfttVB representations (or moments) for Fig. 38 
retrieved size distribution: number-, area-, 

and volume-, log size distributions. 
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4.4 ERROR ANALYSIS OF RETRIEVED AEROSOL 
SIZE DISTRIBUTION 


In order to understand how the (random) errors in optical depth measurunents 
effect the accuracy of retrieved sise distribution results, the ' following 
numerical experiment was performed. Inversions of the simulated T(X) data 
were performed by a nonlinear least squares (NLIS) program. 

In the numerical experiment, the aerosol size distribution n(r) ( cm ^ pm 
was represented by an analytic model, such as, log-normal distribution (LND) with 
three adjustable parameters p^ (i • 1, 2, and 3), defined as (in Ref. 21): 


n(r) 


(37) 


where p^ is the scaling parameter associated with the total number of particles, 

Pj is the mode radius parameter, and p^ the polydispersity parameter. Then, if 
Sto)ces' settling is assumed (Ref. 23), the optical depth T at any time t is defined 


by 


T(t) 


TTL 




r^(t) 


r Q(x,m) n(r) dr 


(38) 


where 


1/2 

r 2 (t) - (o/t) (39) 

Q(x,m) is the Mie efficiency factor, m ■ m' - im" is the complex aerosol 
refractive index, and a is defined in Eq. (51) Section 5.1. 

Values of t were calculated for various times t in the time period 10 
to 3000 secs, for three wavelengths X ■ 0.4416, 0.6328 and 1.15 pm, using 

4 

the following input data: p^ «= 8 x 10 » ?2 “ P3 " 0.5; m ■ 1.47 - 

i(0); path length L » 120 cm; distance from top to laser b^am height, 
h " 100 cm; and particle specific gravity ■ 2.0. The calculated values 
of T(X,t) vs t then provide the simulated (measurement) data. Adding no 
random noise to these simulated data, and using as Inputs in the NLLS code first 
estimates for the three parameters p^ of n(r) different from the values 
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of used to produce the simulated data, we obtain in a few iterations 
the best fit values p.', pi and p* for the parameters. The program also 

A dt <3 

yields the errors p|, p^ and p^ in the retrieved parameters p^. Then 
the error in the retrieved sise distribution n* (r) , obtained by using p|^ 
in Eq. (1), is given by 


An* (r) 
n'(r) 



Ap' ♦ 


A • 

55 - '"“2 


dn* 

¥ 


. Ap* 

3 ■* 


(40) 


The errors in thei* due to errors in p^ are given by the relation 


Al* 

T’ 





(41) 


Next, add random noise to the simulated values of T (t) (obtained by 

4 

using Pj^ ■■ 8 X 10 » ” 2.0, p^ •• 0.5). The random nunber generator is 

used to obtain numbers randomly distributed with a mean value equal to zero, 

and variance equal to unity, i.e., N(0,1). These numbers were then converted 

2 

to a distribution which has a mean y and variance O , by using the following 
relation: 


X - RK *0 + VI (42) 

2 

where X is the random number which is N(V1, O ), and RN is the random number 
which is N(0,1), the number generated by the program. 

It was assumed that the mean of the errors was zero and the standard 
deviation, O, could be expressed as a percentage of the measurement. Thus, 

Tj, - T(l X) (43) 

where is the new measurement with the noise added. 

The values of O used were 0.01, 0.02, 0.05 corresponding to 1%, 2%, and 
5% noise l<*vel. 

Discussion of Results 

The simulated (measurement) data for t (X) were inverted using a non- 
linear least squares (NLLS) program. With no noise added to the simulated 
data, the retrieved size distribution parameters agreed almost perfectly 
with the true values, any discrepancy being negligibly small. 

When noise was added to the simulated data, five sets of noisy data 
were generated for each of the three noise levels (1%, 2%, and 5%), by 
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generating the rAndom numbers five times separately* to allow a bettor 
examination of the accuracy of the inversions of noisy data. Each set 
was then inverted separately. The results are summarised in Table 29. 

For 1% Noise 

For \ ■ 0.4416 bm* the ranges of errors in the retrieved values of p^ 
was 0.66% to 0.84%; in p^* 0.61% to 0.78%; and in p^* 0.74% to 0.95%. 

For X « 0.6328 \m, the ranges of errors in p^ was 0.64T to 0.84%; in 
Pj, 0.55% to 0.78%; and in p^, 0.65% to 0.96%. The range of errors seems to 
widen slightly for the larger wavelengths. 

For 2% Noise 

When the noise in the data increases* T 2 d>le 29 shows that ranges of' 
errors in p^, p^. and for X - 0.4416 pm are 1.30T to 1.66%* 1.21% to 2.11%* 
and 1.46% to 1.92%* respectively; for X > 0.6328 pm* 1.19% to 1.67%* 1.12% to 
1.5S%* and 1.29% to 1.87%* respectively. Again* we notice that for the 
larger wavelengths* the range of errors increases somewhat. 

Errors in Calculated T* vs t 

The retrieved parameters pj were used to plot the retrieved size distributions 
n'(r)* which in turn were used to obtain the calculated (retrieved) T* vs t 
data. Table 30 shows the results of expected percentage errors in the calculated 
values and the ratio of the calculated to the true t values (at selected timen 
in the period 10-3000 sec) for each of the three noise levels. These results 
are for only one wavelength; vis.* X " 0.4416 pm* and for inversion on only 
one set of data. The trends obtained from these results will be typical for 
other cases as well. 

Examination of Table 30 shows that for the three noise levels* the percentage 
errors in the calculated values are generally greater than the percentage 
difference between calculated and true values* indicating that the inversion is 
quite good* even with 5% noise added. One notices that as the noise level 
increases* then both the percentage error in the calculated value and the 
percentage difference between calculated and true values tend to increase. 

The percentage errors in calculated values are least in the time period 
200-1200 sec (i.e.* the r^ region 1.87-4.57 pm). This corresponds to the 
region where percentage errors in retrieved size distributions are also 
minimum and to the region near the peak where the size distribution is changing 
less repidly than for very small size (r -* 0) or very large size (r •♦ “>) 
particles. 
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Errorg in Retrieved Sif Distribution 

From the errors in the retrieved parameters, the errors in the retrieved 
size distribution can easily be calculated by using the Eq. (40) . For data 
with no noise, the errors were negligibly small, as iirould be expected} but 
they become significant when noise is added to the data. Discussions will 
be restricted to the particle size range 0.1 < r < 10 pm. The results are 
plotted in Fig. 41. 

For 1% noise, the error decreases from about 11% at radius 0.1 pm to 1.5% 
at 2.0 pm, then increases to 7.5% at 20 pm. For 2% noise, the error decreases 
from 22% at radius 0.1 pm to 3% at 2 pm, and then increases again to 15% at 20 pm. 
The rate of change of percent error is greater for 2% noise data than for the 1% 
noise data. Thus, as the level of noise in the da... increases, the uncertainty 
in the retrieved parameter P increases and this in turn leads to an increase in 
the uncertainty in the size distributions. 


TABLE 29. RETRIEVED VALUES FOR THE SIZE DISTRIBUTION PARAMETERS AND THE PERCENTAGE ERRORS 
FOR THE NO-NOISE AND NOISY DATA 
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7.879(4) 2.507(3) 3.18 2.06 6.089(-2) 2.96 0.501 1.788(-2) 3.57 99.615 0.017 

7.916(4) 2.918(3) 3.69 2.01 6.749(-2) 3.36 0.475 2.11K-2) 4.44 99.411 0.123 
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TABL£ 30« FIRST SET OF SIMULATED DATA 


X - 0.4416 Mm 


Tine 

sec 

1% 

Noise 

2% 

Noise 

5% 

Noise 

% error 

Calc/True 

% error 

Calc/True 

% error 

Calc/True 

10 

1.71 

1.0012 

3.43 

1.0023 

8.22 

1.0258 

20 

1.69 

1.0011 

3.40 

1.0022 

8.12 

1.0246 

50 

1.5B 

1.0007 

3.18 

l.OOl.* 

7-49 

1.017 

100 

1.35 

0.9997 

2.71 

0.9994 

7.15 

1.0034 

150 

1.17 

0.9952 

2.35 

0.998 

5.51 

0.9941 

200 

1.05 

0.9985 

2.11 

0.997 

4.96 

0.9884 

250 

0.98 

0.9982 

1.97 

0.9964 

4.63 

0.9809 

300 

0.94 

0.9980 

1.89 

0.996 

4.43 

0.984 

400 

0.91 

0.9978 

1.83 

0.9957 

4.27 

0.9846 

500 

0.92 

0.9979 

1.84 

0.9959 

4.27 

0.9879 

600 

0.94 

0.9981 

1.89 

0.9963 

4.35 

0.9930 

700 

0.97 

0.9984 

1.95 

0.9968 

4.47 

0.9988 

SOO 

1.005 

0.9974 

2.01 

0.9974 

4.61 

1.0048 

900 

1.05 

0.9991 . 

2.10 

0.9982 

4.78 

1.0124 

1000 

1.089 

0.9995 

2.18 

0.9965 

4.95 

1.0188 

1200 

1.18 

1.0002 

2.36 

1.0004 

5.31 

1.0323 

1500 

1.35 

1.0133 

2.70 

1.0028 

6.00 

1.0552 

1700 

1.46 

1.0023 

2.93 

1.0044 

6.46 

1.0691 

2000 

1.63 

1.0033 

3.25 

1.0064 

7.09 

1.0875 

3000 

2.29 

0.9954 

4.57 

1.0138 

9.71 

1.1568 

SD 

2.35 X 

-2 

10 

4.71 X 

10-2 

0.106 



Note that in almost all cases scale value agrees with true values to well 
within the error. 
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5 . MODELING OF EFFECTS OF AEROSOL hlCROPHYSICAL AND 
DYNAMICAL PROCESSES ON t ASER BEAM PROPAGATION 


Tha maaBuramants of T(X«t) vb t, shown in Pigs. 3-5 (Saction 2) show that 
tha optical dapth of aarosols variat with tista. Tha raason for thia tiaio bahavior 
of T is that tha aisa distribution n(r) could ba varying in tima. Kara wa assuma 
that during tha coursa of tha axparimant« no chamical raaction takas placa and 
tha changa of tha chamical composition of tha aarosol particlas dua to micro- 
physical procassas is vary small, so that tha rafractiva indax of 
tha particla rnnains constant. Rasults c^tainad frcxs T(X,t) data for diffarant 
timas as dascribad in Saction 4.0 show that the sisa distribution varias with 
tima in such a way that tha size distribution paak valua decraasas with tima 
whataas tha moda radius nhifts toward larger radii. These rasults ara also 
confirmed by tha siza distribution measuraments (Saction 2) . The question 
that arises then is how to explain the tima variation of n(r) and T(X). Tha 
procasses which affist tha time behavior of aerosol particles can be one or 
more of tha following dynamical and microphysical processes i 
Dynamical procassast 

a. Gravitational settling 

b. Stirred settling 

c. Convective Flow due to them:al gradients 

d. Molecular Diffusion 
Microphysical processes t 

a. Thermal Coagulation 

f. Forced Coagulation 

g. Growth/Evaporation 

h. Nucleation 

Xn order to understand how each, or a combination, of these processes 
affect the attenuation of visible and IR laser beams with time, we need to 
perform modeling and sensitivity studies, first for the simple processes, and 
later on, for processes of greater complexity. We decided to model first the 
idealized situations of gravitational settling, thermal coagulation and evaporation 
or growth, for various initial aerosol size distributions. The theories, 
computations and results of these studies are discussed in Sections 5.1 and 5.2. 
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5.1 EFFECTS OF AEROSOL COAGULATION AND SEDIMENTATION 
ON VXSXBLE/XNFRARED LASER BEAMS 


Xntrod^jction 

In order to under etend the eeperate and combined effecta of coagulation 
and aedimentation on the transnittance of laser beams through artificial 
aerosol medium, a joint extinction-coagulation-sedimentation (ECS) model 
has bean set up. This report describes the results of a parametric study 
performed to show the effects of the two processes on the time variation 
of transmittance of visible/XR laser beams through fog oil aerosols settling 
uiK*er gravity in a quiet chmnber. Xn this study, the initial size distribu- 
tions at zero time were varied systematically in such a way that tVieir mode 
diameters lay in the range 0.4 pm to 2.2 pm. Then usino the Mie theory results 
in the ECS model, the corresponding time variation of the optical depth during 
a typical simulated (experiment) run of 0 to 100 min is studied for each of 
the initial size distributions. Preliminary results indicate that for the 
conditions under investigation, coagulation in predcxninant for small particles 
(diameters < 1.4 pm), whereas sedimentation is predominant for large particles 
(diameters > 2.0 pm). 

The attempt was to develop a model that could explain the time behavior 
of extinction measuronents given in Section 2. The theories of differential 
settling and thermal coagulation, described in the aforonentioned references, 
Fuch's boo)( (Ref. 31} and Yue and Deepa)c (Ref. 32) and Deepa)( and Yue (Ref. 33), 
are briefly discussed in later sections. 

Theoretical Considerations 

For the optical geometry of the experiment being simulated in this 
numerical study, consider a box containing a polydisperse aerosol settling 
quietly u,.Jer gravity, through which traverse four visible and infrared 
laser beams, all of them being at the same depth from the’ box top (Fig. 1 in 
Section 2) . 

A plane electromagnetic wave of wavelength X and intensity after 
traversing a distance L through a polydisperse aerosol is attenuated to 
an intensity I, which is given by Bouguer's (Lambert-Beer's) Law, namely 



(44) 


l(X,t) - X^(X) •»pI-T(X,t)] 

O 

wh«r« th« tijiia dapandant optical dapth T(X,t) ia dafinad by 


T(X.t) 


1 

i 


L 

» 

0 




tj'f) t 


-1 t 


0 (>c,ni)iir^n(r»t) dr dl 
*axt 


(45) 


•i 

w)iera t is tha tima> (cm ) ia tha voluma axtinction coafficiantf 

Q (x,m) ia the afflciency factor » x > 2irr/X ia tho aiza paranatar, 

V 

m • m' - m" ia the particle complex refractive index, n(r) ia tha aiza 
-3 -1 

diatribution (cm lim ) , and r^ and are the upper and }ower limita 
of the aeroaol radii. 

Xf we aasune that the aize diatribution remaina unifozn alon? the 
patha of the laaer beona in the aeroaol medium, then 


(X,t) - XL I 

*x. 


r,(t) 


(t) 


0 (x,ro)r n(r,t)dr 


(46) 


We conaider here that changea in the aeroaol aize diatribution and 
hence the optical deptha, with time t occur due to only two microphyaical 
proceaaea, namely, coagulation and aedimentation. Thua the reaulta of a 
joint extinction-coagulation-aedimentation (ECS) model will be described 
here. 

One aim of thia atudy was to determine the conditiona under which 
each of two proceaaea dominates; the other was to develop a auilti spectral 
extinction measurement technique to study the microphysical processes in 
general. For this purpose a numerical parametric study was undertaken which 
involved inputting different starting size distributions n(r,0) at t - 0 
into Eq. (46) and Computing the t (X,t) as f function t for different X. 

Use of analytic model (s) (Ref. 21) for size distribution is best suited 
for such a study. Thus results were obtained with the regularized power 
law (RPL) model which is defined as 
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p(r,0) 


( 47 ) 


Pi Cr/PjT" 


P,-l 


P'5 0^ 

Q (r/p_) ^2 ^ 


For which the mode radius is given by 


■ P2 


f_ 

[l + Pj (p^-1) J 


(48) 


where the parameters p^, p^, and p^ are adjustable constants which can 
be judiciously selected to yield the iuitial size distributions of interest. 

The RPL model is discussed in detail by Deep^k and Box (Ref. 21 ). For r » r^, 
Eq. (47) reduces to a power law, namely, n(r) ^4'*’^^. 

Theory of Differential Settling : Stokes' law states that the terminal 

velocity V of a spherical particle settling under gravity in a quiet medium 

8 

is given by 


2 

9 


% 


( 49 ) 


where p and p are densities of the particle and medium, respectively, 
m p 

r is the particle radius and is the medium viscosity. 

Hence, if a column of uniformly distributed aerosols has a height h 
at t ~ 0, no pu tides of radii greater than r\t) will be present at that 
level any\\rhere in the column at time t, where 

r(t) - (^)** (50) 


and 

« n h 
« 5 ra 

p m 

The upper limit r^ in Eq. (46) is then given by r(t). 
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(51) 




.3 

For alr« tho Bpocific gravity, ■ 1.22 x 10 and viacoslty, 

.422 

■ 1.816 X 10 poiae (dynea -aec/cn ); and g 961 cnt/aac . Thua, for 

-1 2 

oil dropleta, apecific gravity p ■ 0.85, ao thatv^* 1.01925 x 10 r (mn/aoc) 

and for watar dropleta aattling in quiat air, p ■ 1.0, ao that v - 1.1991 
“1 2 

X 10 r (iiun/aec) . Soma valuaa of r and Vg for oil and watar droplata aattling 
in quiet air ara given in Table 31. 

Theory of Thermal Coagulation t When aeroaol particlea come into contact 
and coalesce or adhere to one another, the process is called coagulation. It 
this occurs due to their Brownian motion, then it is referred to as thermal 
coagulation (Ref. 31). A brief theory of Brownian coagulaticxi, baaed on Fuchs' 
treatment, is given aa follows! 

The time evolution of the number concentration N(r,t) having radius r 
in the range r to r ■<- dr at time t, can be expressed as 


TABLE 31. STOKES' VELOCITY FOR SPHERES VARIOUS RADII SETTLING IN QUIET 
AIR IN A CHAMBER. 


Specific gravity ofi 


Viscosity of air, 


air, p^^^ - 1.22 x lo”"' 

fog oil, ■ 0.9218 

-4 2 

■> 1.818 X 10 poise (dynes -sec/an ) 


Droplet 

Radius 

r(ym) 


Stokes ' 

Velocity, V 
5 

Fog Oil Droplets 
(inn/sec) 

Water Droplets 
(mm/sec) 

0.1 \m 

1.019 X 

0 

1 

u* 

1.199 X 10"^ 

0.5 \m 

2.547 X 

io“2 

3.00 X 10”^ 

1.0 ym 

1.019 X 

0 

1 

1.199 X lO"^ 

2.0 ym 

2.038 X 

1 

0 

aH 

2.398 X 10”^ 

5.0 ym 

2.547 


3.00 

10. 0 ym 

10.19 


11.99 
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awtr.t) 

at 


-N(r«t) K(r,r’) N(r*,t) dr* 

^ ^ 1) Ntx^»t) N(r^,t) dr^ dr^ 


(52) 


whcre-X(r^f r^) is the coagulation kernel or frequency of collisions (per 
unit volume) between particles of radiuri r^ and particles of radius r^. 

The first term on the right hand side of Eq. (52) describee the reduction 
of nimber of particles with radius r by the coagulation between particles with 
radius r and other particles. The second term on the right hand side of 
Eq. (52) describes the production of new particles with radius r by the 
coagulation between particles with radius r^ and particles with radius r^, 
so that r^ and r^ are related by the equation 


3 3 3 

IT *f r *■ T 

i j ij 

The coagulation kernel as derived by Fuchs (Ref. 31) is given by 


(53) 


K(ri»r^) - 


°ij^j 


(54) 


where 


’"ij 


- Tf + 


Di + D^ 


"ij 

o 2 ^ 

« (G^ + Gj) 


’j' 






(55) 


and r^ and r. are the radii of particles in the i and j "bins"; 
i 3 

the quantities are diffusion coefficients which can be calculated 

from the Einstein relation, D • kTB, where k is the Boltzmann constant, T, 
the absolute temperature, and B, the mobility given by 
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(56) 


0. &7 


■ “ gjr""- Cl ♦ 1.246 K + 0.42 K • 
6trn r n n 


n 


where n is the viscosity of air, ai..'t K is the Knudsen number, defined by 
m n 

the relation 


i eff 


n 


(57) 


where 


I eff - (Tind^p^^^)"^ 


(58) 


is the effective mean free path of air molecules, n is the number density of 

air moleculesi d ■ 1/2 (d . -f d ,), the d's being diameters of an air 

air aerosol 

molecule and a molecule of the particle under consideration; and, 

”air 

\i ■ ;;; 7-;; , the M's being the molecular weight for air and for 


M , + M , 

air aerosol 

the particle. is the )cinetic velocity of a particle 


./ 


8 KT 

Trm^ 


(59) 


where m^ is the mass of the particle, and 6^^ is a correction factor given by 


3/2 


<60) 


where 


*b 


1TQ: 


(61) 


In order to calculate the change of aerosol size distribution with time 
by computer programming, we have converted the continuous evolution equation 
into an approximate discrete equation by the following formula: 


r. “ r, 2 
1 1 


i (1-1) 


(62) 


where r^ is the radius of the smallest aerosol. Numerical values 2 and -j 
in Eq. (62) are chosen so that the volume of one aerosol in bin i is twice 
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the volume of one aerosol in the preceding bin. In our model, we have used 
35 bins to cover the aerosol sise spectrum from 0.01 pm to about 25 pm. 


Computations and Results 

A parametric study of the separate and combined effects of sedimentation 
and coagulation on the optical depths for four wavelengths (X - 0.4416, 
0.6328, 1.15, and 3.39 pm) has been carried out for fog oil aerosols, using 
the RPL size distribution model in Eq. (47). The input data for the 
physical properties of fog oil, including the real and imaginairy parts of 
complex refractive index for fog oil at the four wavelengthtt, are given 
in Table 1, Section 2; h« 62.5 cm and L ■> 117 cm. Optical extinction 
values were computed with the use of Mie theory cMnputor codes. The geometry 
of the simulated optical system is shown by the schematic illustration in 
Pig. 42. 


In the model, we assume different initial mode diameters 


d at t 
in 


0 by choosing parameters p. 


p^ " 4.0 and different values of 


the parameters p^. Parameter p^ is the scaling parameter which is obtained 
by normalizing T ■ 7.0 for X “ 0.6326 pm at t 0. The choice of the 
value 7.0 for normalization purposes has no special meaning other than the 
feet tl.at it is- representative of the experimental conditions under investi- 
gation. Computations were carried out for different initial size distri- 
butions with a different mode radii obtained by choosing different sets of 
p^. For clarity, the results of only three sets of p^ values (Table 32) are 
discussed here. The size distributions obtained at times t • 0, 1.5, 25, 50, 

75 and 100 min during a simulated (experiment) run are graphically illustrated 
in Figs. 42 to 46. Figures 43, 45, and 47 are for coagulation and sedimentation 
combined and Figs. 44, 46 and 48 are for coagulation alone. For sedimentation 
alone, the shape of the size distribution should remain the same as the 
initial one at t ■ 0 min, except that the upper di^uneter "cut off" (d 

max 

continues to decrease as time increases. For example, at the times given, 

the d » 20.48, 15.96, 3.91, 2.76, 2.26, and 1.95 pm, respectively. These 
max 

values correspond to a beam (Fig. 42) at h « 63.5 cm and fog oil droplets with 

-4 

p ■ 0.9214 settling under gravity in air with viscosity » 1.813 x 10 
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FIGURE 43. Aerosol size distributions st times t " 0, 1.5, 25, SO, 75, 100 min due 
coagulation * sedim^tation. Initial mode diamter d • 1.6 ]m. 
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FIGURE 44. Aerosol size distributions at times t 0, l.S, 25, SO, 75, 100 min due 
coagulation alone. Initial mode diameter d • 1,8 pa. 




FIGURE 45. Aerosol size distribations at tines t ■ 3, 1.5, 25, 50, 75, 100 min due to 
coagulation * sediments tim. Initial aode diaamter d • 1.0 ys. 
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TABLE 32. SUmtARY OF THE PERTIHENT VARIABLES INVOLVED IN THE PARAMETRIC NODELING STDDY 
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*”** Pp, ■ 1*225 X 10 (Bq. 49). Th« results lor the time dependence of T(X), 
obtained for the different initial sise distributions by the use of the ECS 
model are shown in Pigs 49 to 57. The results for the combined coagulation * 
sedimentatiOT effects are illustrated in Pigs. 49, 52, and 55i while those due 
to coagulation and sedimentation separately are shown in Pigs. 50, 53 and 56, 
and 51, 54, and 57, respectively. 

Discussion of fiesults and Conclusions 

Prom this study it becomes evident that using the RPL sise distribution model, 

the relative strengths of the coagulation and sedimentation effects on T vs t 

behavior depend on the initial mode diameter d^^ such that coagulation effects 

dominate for d less than about 1.0 to 1.4 ym; while sedimentation effects 
m 

dominate for d greater than about 1 .6 ym* 

Rl 

The initial aerosol si«e distribution with the largest mode diameter 

(vis., d||| ■ 1.8 ym) is shown in Pig. 2. A comparison between figs. SO and 51 

shows that the sedimentation effect is dominant, since the coagulation effect 

causes a slow time variation of T for all four laser beams Cfig. SO). But, 

when the initial d ■> 1.0 ym coagulation effects on T(X,t) neem to become 

m 

comparable to sedimentation effects, as shown in Figs. S2, S3, and 54. 

Pigure S2 for the combined effects shows that for X ■ 3.39 ym, T initially 

increases for the first 20.5 min, and decreases thereafteri for other 

wavelengths, it simply decreases with time. Por coagulation alone T pea)cs 

at 49.5 min for X ■ 3.39 ym (Pig. 53). 

When d *0.4 ym, the sedimentation effects on T vs. t are insignificant 
m 

(Figs. 55, 56, and 57) i but due to the increase in small particles, coagula- 
tion beccxnes increasingly dominant, thereby causing the time rate of change 
of d to increase. This in turn causes the time variation of T to differ 

PI 

for the differing wavelengths! w.g., for initial d^ - 1.0 ym, the T vs. t 

curves depict a modal behavior for X ■ 3.39 ym, but not for the other three 

wavelengths. But as the initial d^ is reduced, the modal behavior of T vs- 

t curve becomes extended to smaller and smaller wavelengths. Thus, when 

initial d ■ 0.4 ym, T vs. *- curve for X ■ 1.15 ym dejpicts a model behavior 
m 

as well. 


120 




FIGURE 49. Optical depth T vs. time t data for wavelengths 0.4416, 0.6328, I.IS, 3.39 Vm for 
coagulation + sedimentation. Initial mode diameter d * 1.8 Ub>* 
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FIGURE 52. Optical depth T vs. time t data for wavelengths 0.4416, 0.6328, 1.15, 3.39 ya for 
coagulation + sedimentation. Initial mode diameter d • 1.0 lua. 
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FIGURE 57. Optical depth I vs. tine t data for tavelenpths 0.4416, 0.6328, 1.15, 3.39 Un for 
sedimentation alone. Initial mode diameter d ■ 0.4 yn. 



The results described here arc only preliminary and are based on only 

one sise distribution model and in which only the mode diameter was varied. 

Even on the basis of these somewhat scant results, one can easily conclude 

the relative importance of coagulation and sedimentation mechanisms depe lus 

on several parameters, such as, wavelength, X, and h (Pig. 42), end relative 

concentration of smaller size particles given by the mode diameter d_. 

in 

Such a parametric study, in addition to providing an understanding of 
the effect of coagulation on laser beam propagation in aerosols, can lead 
CO the development of simple and accurate experimental techniques for 
investigating the average microphysical processes. An important conclusion 
that can be drawn feexn such a study is that by performing the extinction 
experiment for a single wavelength in the visible region, one can easily 
miss noticing the effects of processes like coagulation. Therefore, it is 
recommended, as was proposed in Ref. 32; that in order to get a more realistic 
picture, one must perform these measurements at multiwavclengths, preferably 
spanning a wide spectral range, from the visible to the IR. Work is under 
way to extend this parametric study to include stirred settling effects of 
polydispersity and number concentration, particles of different chemical 
composition (refractive index) and shapes. 
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5.2 MODELING THE TRANSMISSION OP VISIBLE/IR LASER 
BEAMS IN POLYDISPERSE AEROSOL PARTICLES UNDERGOING 
COAGULATION AND EVAPORATION OR GROVH'H PROCESSES 

According to the wnll-known Bouguer's (Lainbort-Baor'n) law, tho ratio of 
the intenaitiaa of a plane electromagnetic wave after and before traveraing a 
distance through a polydisperse aerosol is a function of tho optical depth of 
the medium. The value of this optical depth depends on the size distribution 
and optical properties of the aerosol particles under consideration. If we 
assume that during the course of the experiment, no chemical reaction takes 
place and the change of the chemical composition of the aerosol particles due to 
microphysical processes is very small, we can assume that the refractive index 
of the particle remains constant. However, the aerosol size distribution will 
be inadvertantly changed and consequently the transmission of visible/lR 
laser beams in polydisperse aerosol particles will show a time-dependent 
behavior. The processes which will affect the size distribution of aerosol 
particles can be divided into two categories: dynamical processes and micro- 

physical processes. The name of each process and condition, under which size 
distribution will be significantly affected are outlined as follows: 

Dynamical Processes 

a. gravitational settling — will affect the larger particles 

b. stirred settling- -particles of all sizes will be affected depending 

on the time duration and strength of stirring 

c. free convection due to thermal gradient — particles of all sizes 

will be affected depending on the magnitude of the thermal 
gradient 

d. molecular diffusion — depends on the total area of the enclosure 

exposed to the aerosol peurticle. 

Microphysical Processes 

a. thermal coagulation — depends on three factors: time, concentration 

of particle and spread of the particle size distribution 
The longer the experimental time, the higher the concentration 
of particles and the more spread of particle size i^pectrum will 
m^d(e the coagulation process more dominant in changing the 
aerosol size distribution. 
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b. forcad coagulation — dapands on whathar axtarnal forca othar than 

gravitation ia praaant 

c. growth or avaporation — dapanda on tha physical propartias 

(tang>aratura* vapor praaaura) of tha particla and air parcala 
surrounding tha particla. 

d. nucleation — dapands on tha concantration of gasaous spacias in 

tha ambiant (undar ordinary conditional this process can bo 
ignored. ) 

Sinco mora than ono procass will ba taking placo at tha sama tima and diffaront 
procassas will bo coupling with aach othar, numarical simulation of tha changa 
of optical depth for aerosols undergoing different dynamical and microphysical 
procesaes at the same t.'ma will ba very complicated and the affect of each 
procesa to the tremsmission of laser beams can hardly be interpreted. The effect 
of each process will be studied separately and in this section, we concentrate on 
th'i study of the effect of evaporation/growth on the temporal changa of optical 
depth. Since the data for vapor pressures of fog oil were not available, we 
proceeded to obtain numerical results for water fog aerosols, which contain 
solution of sea salt. 

We assume that tha amount of sea salt particles in each particle is such 
that tha water vapor pressure above the surface of the particle is in equilibrium 
with that in the anbient. We have the familiar expression 


S 



where 


A 


20 ' 

P'RT 


( 63 ) 


( 64 ) 


B 


i p. M 
d w 

M p 
s w 


( 66 ) 


where S - relative humidity, r - radius of droplet, r^ - radius of dry salt 
O' «■ surface tension of the aqueous solution 
p' ■ density of the aqueous solution 
R - universal gas constant 
T <■ absolute temperature 
i “ Van Hoff factor 
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■ danslty of th« dry particles 
- density of water 

■ molecular weight of water 

" molecular weight of the soluble salt 

The amount of sea-salt in each droplet can easily be calculated by 
applying Eq. (63). As soon ai> the relative humidity of the ambient air is 
changed, the equilibrium radius will be changed accordingly. In order to 
• estimate the effect of relative humidity on the size distribution of fog, we 
have used the measured size distribution of fog at R.H. 99.5% to calculate 
the value of r^ for each value of r at that relative humidity. We assume 
that the dry salts are sodium chloride particles. Then we calculate the size 
distribution at an environment of R.H. ■> 98.5% and 97.5%. The measured size 
distribution at 99.5% and the calculated size distributions at R.H. ■ 98.5% and 
97.5% are shown in Fig. 58. It can be seen that a small change of relative 
humidity will cause a significant change in the aerosol size distribution, 
especially for those aerosols with larger radii. 

In our model studies, we assume that the aerosol particles in the enclosure 
are sea-salt containing particles whose initial size distribution can be 
described by the regularized power law (RPL) model, defined by Eqs. (47 and 48) 
in Secticm 5.1. 

We assume p^ « 4.0, p^ - 4.0 and the mode radius r^ - 1.4 pm. From 
Eq. (48) , the value of p^ can be calculated. The value of p^^ is adjusted 
so that tiie optical depth at X ■ 0.6328 pm is 7.0. We arbitrarily let the 
relative humidity of the ambient air change at a rate of 1% per hour. The change 

of the optical depth with time for different wavelengths is shown in Figs. 59(a) and (b) 

From our preliminary results, it can be seen that even a very slow change of 

temperature will cause a great change of optical depth. This is primarily 

due to the fact that equilibrium sizes of solution droplets are veiy' sensitive 

to the relative humidity of the ambient air. It is expected that if the 

aerosol- particles under consideration are oil fogs, the effect of condensation/ .. 

evaporation on the optical depth will not be so significant. 

The computer programs to study the combined effects of coagulation and 
evaporation on laser beam propagation have been developed, but have not been 
optimized as yet. Work is still in progress. Therefore, the results for 
the combined effects of coagulation and evaporation are not shown here. 
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FIGURE 58. Measured (at R.M. * 99.5%) and calculated (at R.H. - 99.5 
and 97.5%) size distributions. 
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Figure 59(b). Optical depth versus time for particles undergoing 
evaporation when the relative humidity changes by one percent per hour. 
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6. MULTIPLE SCATTERING EFFECTS OF LASER BEAMS 
TRAVERSING DENSE AEROSOLS 


In recent years there has been a considerable increase in interest in the 
problem of multiple scattering effects on laser beam propagation in an aerosol 
media. In two recent papers (Ref. 34, 35), we have investigated the case of 
forward (single) scattering as it effects the transmittance or extinction 
coefficient measurements for plane parallel (collimated) beam of radiation 
traversing aerosol media. In an earlier work (Refs. 5 , 6), we used the 
successive scattering approach to determine the beam broadening effects of 
second and third order scattering on a collimated pencil of radiation (white 
light) . Since then several papers have appeared in open literature on the 
subject of aerosol effects on laser beams. 

Most of the approaches adopted in MS studies of beams of radiation 
are based on the assumption of small-angle approximation. They can be categorized 
as either analytic approaches (Refs. 36-46) or successive scattering computer 
simulation approaches (Refs. 47-51). Some of these studies deal with incomplete 
MS problems, namely, with two or three orders of scattering (Refs. 6, 52-55. 

Among the analytic methods, most of the approaches are based on the 
solution of the radiative transfer equation (RTE) for narrow collimated beeuns 
traversing dense particulate (aerosols or colloid) media (Refs. 38, 41, 56-61). 

An exact solution of the RTE for beams traversing media with highly anisotropic 
phase functions, is possible only by making the small angle approximation. 

However, to obtain quatititative results, one needs to make further approxima- 
tions, each leading to markedly different results. A brief review of three 
approximations will be discussed in the next section. 

If a relatively narrow beam propagates in a scattering medium, photons 
are constantly removed from the beam. However, if the scatterers are of a 
size equal to or greater than the radiation wavelength, such as in the 
case of smoke, dust, or fog particles compared to visible wavelengths, then 
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most scattaring sventB will raault in a ccHigiaratively amall deflection of 
the photon. Thia may lead to a gradual apreading of the original beam, both 
in thickneaa and angle. 

Zn thia report, we exan..ne the aignal that may be detected, aa functicma 
of both experimental geometry and the propertiea of the acattering medim. 

We ahall employ the small -angle approximation to the equation of radiative 
transfer, which ignores photons which have suffered large deflections as 
they will be assumed lost. In order to obtain tractable answers, it will 
prove necessary to assume simple analytic forms for the scattering phase 
function and the initial beam profile. Nevertheless, the analysis presented 
in this report will be as free as possible of unnecessary approximations. The 
work discussed in Section 6.1 has been reported in Ref. 45. 



4 

li 

6.1 GENERAL SOLUTION OF RTE IN THE SMALL ANGLE APPROXIMATION 

Let I(Zf 0) dV dfi be the intensity of radiation (or the nunber of 
photons) in a volume element dV centered at the point z« r ■ (x« y) » and 
travelling within a cone of solid angle dn centered about the direction 0. 

Then I satisfies the radiative transfer equation, which we may write 
** (Refs. 43 and 44). 

fi.^ + 0 l(z,r, n)-u>o 

where o is the extinction coefficient (km '*') 
is the albedo of single scattering 
\l> = cos ^ (n . u') is the scattering angle 
and P(lp) is the scattering phase function (sr ^) 

Even allowing for cylindrical symmetry about the axis of propagation 
(the z axis) , Eg. (66) is exceedingly hard to solve numerically. However 
since the diameter of our detector will always be small compared to the total 
propagation distance, we may safely assume that all photons which are eventually 
detected will have spent their flight time travelling essentially parallel to 
the z-axis. We may thus set cos 6 1, where 6 is the angle the photon makes 

with the z axis. 

Note that this approximation ignores the contribution from all photons 
which undergo at least one large-angle scattering event. All such 
photons will clearly need to undergo at least a second large-angle scattering 
event, and maybe even a third, in order for them to reach the detector. As 
we are assuming that the phase function, P, is strongly forward-peaked, the 
probability of two or more large-angle scatterings is clearly very small, 
and thus the neglected contribution will be small. 
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The main effect of this assumption is to replace the unit propagation 


vector by 


n = (n , n )■► { 1 , n ) 
.. z -X -X 


(67) 


Although this new propagation vector is no longer correctly normalized, this 

should not cause any problems, as the number of photons for which |n | « 1 

-.L 

does not hold, will clearly be small. 

The second effect is that we may use n, - n' as the argument of P in 

-X *1 

Eg. 1, i.e.. 


P(n«n')-> P(n - n') 

~X 


(60) 


and third effect is to replace the limits of this (two dimensional) integral 
by + ®. With these points in mind, we may now rewrite Eq.(66) as 


3 3 f f ^ 

(3? * £' - “o " J ■ 2i’ si’ '*=i 


Equation (69) is referred to as the radiative transfer equation (RTE) in the 
small-angle approximation. Its main advantage over Eq. (66) is in the simplification 
of the directional derivative. This equation has been used extensively in the 
theory of foil penetration by fast charged particles (Refs. 62-64) . Though 
Wentzel (Ref. 3) was the first to use the small-angle approximation for charged 
particle transfer, the first person to employ this equation in the field 
of radiative transfer appears to be Dolin (Ref. 59). 

One further result of the small-angle approximation is that all detected 
photons are assumed to have travelled the same distance. Thus their time of 
travel is constant, and 'a pulse will undergo no time -dispersion. 
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Fornml Solution 


Equation (69) may be solved, at least formally, by the use of Fourier 
transform techniques. Introducing the definitions 

f(z, n. C) - (271)'^ |*||j I(z, r, n ) e^^D’£ * dr dn^ 


(70) 


and P(C) » (211)’’ I I P(n ) o ~ '-L dn, 

' } 1 -1 •'1 
•OD 

we tfdce the double Fourier transform of Eq. (69) to obtain (Ref. 59) 


(71) 


[h " D . ^ + oj i(2» n» U - 211 u)^o p (O i(z, n,C) 


(72) 


Equation (72) is easily solved, to yield (Refs. 59, 64). 


1 (z, n>C) 


1 (n, C + * 

O ^ 


n) e 


- oz + n 


(73) 


where 


n = n(z. n,U 


2n U) o 
o 



+ z' 


n I ) dz • 


(74) 


A 

and lj^(n» C) is the Fourier transform of the initial intensity distribution 
(incident beam profile) at z <■ 0 . 

To obtain the intensity distribution at any point in the medium, it is 
merely necessary to re-transform Eq. { 73 ) (Refs. 64 emd 65), 
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-2 f f f f A - i(n . r ♦ C . n ) 

Kb, r, n^) - (2u) J J J J ^<*» !}» §) • dn dC (75) 

where ^ is given by Eq. (73 )• 

In Eqs. (74) and (75), there are five integrations to be performed in 
order to obtain numerical values of the specific intensity I(b, r, n,). In 
principle, it should be possible to perform these integrations numerically; 
but in practice it is prohibitively expensive to do this using computers 
such as CDC 6600 and 7600, due to the necessity of evaluating the 4- 
dimensional Fourier transform in Eq. (75). 

One way to make the problem tractable is to convert the integration over 
the rapidly oscillating function e»qj(-i(n • r -f C * n,)J in the Fourier trans- 
form into integration over a distance x of a slowly varying function, even 
though multi-dimensional. Such an approach was adopted by Tam and Zardecki 
(Ref. 43) for the case of a Gaussian phase function. They expanded the expo- 
nential facto., containing 0 into a power series, in which each term could be 
analytically Integrated over the four variables (n, C) • The multi-dimensional 
integration over z then has to be done numerically by employing any of of the 
standard procedures, such as, method of Lyness (Ref. 66), Conroy (Ref. 67, 68), 
Monte Carlo (Ref. 69), etc. The details of the method are described in a later 
section. It should be noted that the approach of Tam and Zardecki has not yet 
been tried for non-Gaussian phase functions. 

Another approach, which is approximate but considerably simpler since 
it avoids multi -dimensional integration, was used by Arnush (Ref. 38) and 
Stotts (Ref. 61) in solving the problem of a radiative transfer of a 
collimated beam. Arnush (Ref. 38) treated the case of a narrow collimated 
beam in ocean, and Stotts (Ref. 61), that of a laser beam traversing 
particulate medium. In both cases, they solved a non- homogeneous RTE (con- 
taining a source term) under homogeneous boundary conditions, in this report, 
we consider a homogeneous I^E (in which source term is absent) and solve 
it under non-homogeneous boundary conditions. As is well known (Ref. 70) , 
the two formulations are essentially equivalent. Therefore, their method, 
which will be referred to hereafter as the Arnush-Stotts type (AST) 
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«csa.«s.“! 



approximation, could be adapted to our case. 

The AST approximation basically involves a series expansion of the 

integrand, in R (Eq. 74), and retaining the first two terns, which subsequently 

enables one to perform analytically the 4-dimensional integration in Eq. (75). 

The details will be discussed in a later section. 

There is yet another approximate method, which avoids the difficulties of 

carrying out multi-dimensional integration numerically. This method, hereafter 

referred to as the Dolin-Pante method, essentially involves the expansion of 

I(z, r, n'_^) in the integrand of Eq. (69) in a Taylor series about n^ ■ nj^. 

Subsequently, the RTB in the small-angle approximation (Eq. 75 ) is replaced by 

a system of two equations — one for the unscattered radiance 1^°^, and the 

other for and the scattered radiance Once the equation for 

(s) 

is solved, I can be found by solving a non- homogeneous RTE, with the 
source term being determined in terms of 1^°^. The details of the method are 
explained in a later section. 

In the work reported here, we evaluated the detected power P for a 
6-function laser beam propagating through dense aerosols for different 
phase function models. The other quantities of interest are irradiance 
and the beam spread. For which the formulations are given below but no 
computations have been performed as yet. 


Quantities of Interest; Irradiance , 

Power Detected and Beam Spread 

Irradiance : In this study, if is the irradiance (flux density) and 

received power, rather than radiance, which is of interest to us. Using the 
relation between irradiance, N, and radiance, I, we may simplify Eq. (75) 
(Refs. 59 and 65). 

N(z, r) = I l(z, r, n )ln . z] dn 


I(z» r, 




n ) 



(76) 



n. 0) 


e ~ ~ dn 


(77) 
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With th« •limination of C in Eq. (77), wo nay aiuqiltfy Eqs. (73) and (74) 


i(x, n, 0) 


^ -OB ♦ 

I (n, B n) • 

o - - 


(73*) 


whara 0 

o 



n|) dB* 


(74*) 


Equation (77) can now ba furthar ain?>lifiad by an appeal to ayininetry. Sinca 
P(ij») clearly depends only on the scalar | r^ | , and not on the vector P 

P will similarly ba a scalar function, as will 0^. Similarly, if wa assuna 
that the incident beam profile is circularly symmetric, then lQ(n, b, n) will 
also be a scalar function of Thus X(b, n, 0) will ba a scalar function, 
and Eq. (77) becomes 


N(z, r) ■ 2n 


m 
• n 


Jo(n 


r) Kb, n, 0) n dn 


(77* ) 


From Eq. (77*) we see immediately that N is a scalar function of r, as we 
would expect from the above symmetry arguments. A more ‘ tractable expression 
for the case of 6-fun<:tion will be given later in Eq. (96) . 


Detected Power t In most instances, of course, what we are most interested 
in (and what we physically measure) is the power received by some detector. 

This will involve the integration of Eq. (77*) over the area of the detector, 
perhaps modulated by a response function. If we assume a coaxial, circular 
detector of radius R, with a flat response, then we have 


P(z, R) - 2lf 



r) 


r dr 


4Tt 


"Oz 


jj^(n R) 


io(n. 


z n) e R dn 


(78) 


This result should prove' amen<able to numerical integration, especially if a 
relatively simple expression for (2^ can be obtaii ed. Seunple results will be 
presented. 
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One useful result which can be obtained analytically, is the total power 
crossing a surface z " constant! 


P{z, “ I I V **£ 


- I |*f I 2^ e"^ 5 ' £ dn dr 


■ 4n I(s, 0, 0) 


_ ^ -oz + 2v u) a z P(0) 

471 i (0, 0) e ® 

o 


■ F e 
o 


-(1 “ u) ) oz 
o 


(79) 


where F is the incident total power, and we have used the fact that P(0) ■ (2 tt) 
o 


From Eq. (79) we see that the only energy removed from the beam is that lost 
by absorption, i.e,, there is no bac)cscatter , 


fleam Spread : One further parameter which will often prove useful is the 

beam spread, which we may define as 


<r^ 


r) r dr 


£ N(z, r) r^ dr / fo N(z, 

r N(z, 

Jo 


271 F 


-1 


(1 > ( 1 ) ) oz 
o 


r) r dr 


(80) 
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EXACT SOLUTIONS OF RTE IN SMALL ANGLE APPROXIMATION 
FOR GAUSSIAN BEAN AM) VARIOUS PHSE FUNCTIONS 

Xn this ■•ctlon, w« shall dsvisa axpraaaiona for tho two quantitloo” 
namaly* P and <r > — fron the solution of ths RTE in tho snail angla mpptoxi^ 
nation, for ths caso of a lasor beam, with a Gaussian intonsity profile, 
traversing a scattering medium, represented by four simple phase function 
models— *one Gaussian and three non>Gaussian (vis., two exponential and 
binomial) functions. An alternative expression for the beam profile N 
will also bo given. 

The power detected will be obtained for a collimated beam with zero 
width in space, and the expressions for beam spread, for any realistic beam. 

At the entrance (z ■ 0) to the scattering medium, a laser beam profile is 
represented by a Gaussian functional form, both for the radial intensity 
distribution, and the angular divergence, given by 

I_(r, n ) - F e)g>(-B^ nf - Y^ r^) (81) 

o ~ -j. o 'J. 

This may be easily transformed, and, in particular, we have 

i_(n» z n) * ^^(2w) ^ exp(-n^ / 4 y^ " nV 4 b^) (82) 

o o 

Xn general, the laser beam profile will be well-collimated, so that B and 
Y will be large. The inclusion of Eq. (82 ) in Bq. (78) will in no way com- 
plicate the numerical integration, though in our examples later we will allow 
both to go to infinity, so as to reduce the nvu.iber of parameters whose 
Influence should be examined. Xn all practical calculations, however, realistic 
values of both parameters should be included. Next, we obtain the three 
quantities for the four phase function models. 

Gaussian Phase Function 

A Gaussian functional form is also often en^loyed to describe P(t)^), since 
it can approximate the forward peak of the scattering pattern for a spherical 
particle on which a plane wave is incident. Thus, 


wh«r« a is an adjustabls parameter « which controls the shape of the forward 
peak and is related to the nss scatterin 9 angle )|) defined as 


• 2V 


I P0j>) d<» 


(83a) 


2 -2 

Jt is easily shown that for the Gaussian case ){> • a . (Though a will 

usually be large, it will rarely, if ever, be as large as 6 or 

Taking the Fourier transform of Eq. (83), we find 


- I J (C H/) P(l(») ^ ^ 
in O 


(71*) 


-CV4 


/ 211 


(84) 


and hence 0 


0) o n erf (zn/2a) 

o 


(85) 


Substituting Eqs. (16) and (19) into (12) we find 


p(r, R) - P R I J (n R) ejq> [u) 0 n ^ o erf (zn/2d) 
o Jo 1 o 


- oz - n^/4 n^/4 dn 


( 86 ) 


Equation (86) is an exact equation, which can be solved numerically for 
various values of 3 and y. In order to reduce the number of parameters, we make 
a simplification in Eq. (86) and consider the limiting case of B, Y This 

physically implies that the beam is collimated euid has zero width in space. 

Then, by making the variable changes 

T - 0 z 

T » 0) T 

s o 

2 1/2 

G - R/z 

X - HR 
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Equation (86) can bo raducad to 


Pit, R) - r 

o 


■’f 

# A 


<X) *xp /n Gx"^ arf (x/2C)] dx 


(87) 


T)ia powar of t)ia unacattarad baan at an optical depth of T ia, of couraa« 
P^e ^ . Thua, the praaanca of forward acattaring haa increaaad the datactad 
power by the amplification factor A# given by 

A(T^, G) - C J^(X) exp jr^ /n G x"^ erf(x/2G^ dx (B8) 

where A is a function of two parametera; the acattering optical thic)inesaf 

and G, the geometry factor. 

Finally, we may obtain the beam spread for the general case of a beam 
given by Eq. (81) by substituting from Eqa. (77) and (82) into Eq. (80): 

<r^> - i z^/3a^ + 

■ (89) 

In general, the first tenn should dominate, except perhaps, close to the 
point of entry into the medium. 

Non-Gausslan Phase Functions 

Although the Gaussian form in Eq. (83) is a popular model for the for- 
ward pea)( of the phase function, it is often a good idea to examine other 
models, to ma)ie sure that none of the results are singly an artifact of the 
Gaussian model. In this section, therefore, we shall examine a number of 
other functional forms which may be (and have also been) used to model aniso- 
tropic phase functions. We shall follow essentially the same steps as in the 
previous section, and present only the results, unless further explanation 
is necessary. 

Exponential Phase Functions : 

For 8 infinitely narrow beam 

2 -Ct ill 

P(l|i) - a e /2H (90a) 


148 




/ 2TI 


(90b) 


P(C) - a^(a^ ♦ 
n - T. (1 ♦ 


where 


y ■ z r\/a m X/ G 


Thus 


^QD MM 

A(T^, G) - j J^(X) exp jr^d + xV6<5^)"^^J 


dX 


(90c) 

(91) 

(90d) 


and for a general beain 

<r^> - 2 T zVa^ + /&^ •*■ y ^ 

0 


(90e) 


Phase Function for Sea Water : Similarly, for the infinitely narrow 

beam 


P({p) m a '^/2 ti 


(92a) 

P(0 - a(a^ + 2ir 


(92b) 

fio ■ Tg y"^ In jy/*^ + (1 + 


(92c) 

A(Tg,G) - Jj^(X) 

In [x/G /2 + (1 + 1 xV^)^''^J dx 

(92d) 


and, for a general beeun 

<r^> - 2 T z^ / 3 + z^ / 6^ + Y’^ (92e) 

8 


Note that, although Eq. (92a) implies P(0) ■ «*», 
solid angle factor, leads to a finite result for 
through any angle. In fact, Eq. (92a) has been 
et al. (Ref. 65) to model the phase function of 


the inclusion of the correct 
the amount of light scattered 
employed by Bravo-Zhivotovs)ciy et al . 
sea water. 
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Binomial Phase Functions: 


This time, wo consider phase functions based on the functional form 
(1 4 We will need the result that 

aCD 

I Jp(nil') (1 ♦ d|f» - (n/ 2 o)^ K (n/a) / 4 1 ) (93) 

)o 

where is the modified Bessel function of the second kind. Thus if 


P^(4.) 


2 VI a U 


4 aV)"^'^ / 


2n 


(94a) 


f^(U - (C/ 2 a)^ (C/a) r (H) 

flo - ig /n r(y 4 i) [k^ (yM ^^.^(y') + K^_i(y) i'y(y')] / T(vi) 


(94b) 


(94c) 


where y* = y/p-l 

and is the modified Struve function of order Vi« 

The expression for A may be easily written down. For VI > 1» we may 
obtain the beam spread for a general beam ($ , y finite) 

<r^> » T z^/3 a^lu - 1) 4 + y"^ (9^^ 

s 


Note that if V* is .an odd half-integer , Eq. (94c) may be expressed in terms 
of exponential functions. For example, for p » 3/2 we find that for 
infinitely narrow beam (6, r ■♦■ “) 


n (3/2) - T jj /2 y~^ - (14 2/2 y"^)J (95) 

o s 

A(Tg, G) - r Jj^(X) expjx^ [2 /2 Gx“^ - (14 2/2 Gx"^)]j dx (95* 
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Beam Profile U(r) 


The irradiance profile of the expanding beaip is also a quantity of interest, 
even though no numerical results are presented here. The formal expression for 
the transverse irradiance profile is given by Eq. (77*). In most cases, this 
expression is well behaved. However, for the special case of 6» Y Eq. (77') 

will diverge. Therefore, an alternative expression for N(r) may be obtained 
either by integration by parts of Eq. (77') or by differentiating Eq. (78). 
Adopting the second approach and setting I « P (2ii) , we obtain 


N(r) 


1 ^ 
2TIR dR 


1 T / ‘ 

= - F^(2n)" r" e' J^^(X) e ° y dX (96) 

o o 

where the prime denotes differentiation of with respect to y (Eq. 91). 

With the exception of the sea-water phase function, (2^ goes as y ^ for large 
y, and so (2* goes as y and convergence is assured. 


THREE ALTERNATIVE APPROACHES FOR SOLVING THE RTE 

For the sake of completeness, we shall give the mathematical formulations 
for three methods — one exact (Tam and Zardecki method) and two approximate 
(Arnush-Stotts type and Dolin-Fante) methods. However, computations will 
be performed for only the Arnush-Stotts. 

Tam and Zardecki *s Exact Method 

The method ot xam and Zardecki is exact, at least in principle, but 
requires the evaluation of multidimensional integration, the order of which 
is equal to the order of multiple scattering required. We will restrict this 
discussion to the case of the Gaussian phase function only. 
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The Tam and Zardecki method consists in cxpandino exp ({] ) in a Tay*or 

o ^ 

series, before performing the integration over z' <Eq. 74'). Thus inserting 
Eq. (84) in Eq. (74*), and performing the Taylor expansion, yields 


exp 


(fl ) - 
o 


1 + 


I 


ma 


^ f-'-r 

^ i*" m! ^ ^ 


dz. 


dz^ exp { 


m 


m 

i£i 


zl / 40*) 


(97) 


Wa may now perform the inverse Fourier transform (Eq. 77') to yield 


m 


N(z, r) * F e”^ (Ti z*)"^ “fr N (z, r) 

o m=*0 m ; m 


(98) 


2 «2 2 2 .2 ^2 
where N = exp {- 

zy+B zy+e 


(99a) 


and 


N = [ ■ * ' f dZi * ' ‘ dz A exp f-r^/z* A 1 

n> Jq Jo ^ ^ ^ mi 


(99b) 


where A 


m 


-2 ”* 2 - 2-2 -2 


(99c) 


We may note in particular that may be evaluated analytically in terms 
of the error function. The resulting expression is quite complicated, except 
in the case where 3 and y go to infinity, in which case we get 


where 


Nj^(z, r) = a* /ir [i - erf(g)] / 2g 
g * r o/z 


(100a) 

(100b) 
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Turning our attention to the power received# we obtain an expression 
similar to Eqs. (98) and (99) , viz. 


.m 


P{z, R) 


-T 


''o • ■ -fr 'm'*- 


where P 


and 


«2 2 

1 - exp {- — — - — ^—) 

Y + 3 


m 


f ••• f dz ••• dz^ {l - exp [-rVz^ A ]} 


(101) 

(102a) 

(102b) 


Note that, from Eqs. (101) and (102), Eq. (79) may be obtained trivially. 

As with P^ is also analytic , and in the simple case of B, Y "*■ “» obtain 

^2 

P^U, R) » 1 - e + G /n [l - erf(G)J (103) 

The number of terms required for the convergence of the ser'ies in Eqs. (98) 
and (101) grows steadily with T , and so in some case for large optical thicknesses 
it may become prohibitively expensive to use it. Nevertheless, these results 
have one use in that Tam and Zardecki (Ref. 43). have shown that the mth order 
terms in Eqs. (98) and (lol) correspond to the contribution from mth order 
scattering. This in itself is a useful result. 

Another use suggests itself, however. The Gaussian phase function is simply 
a model, with the parameters ot and ii)^ available for adjustment to match real 
scattering patterns. Since we now have a simple expression for the singly - 
scattered contribution, we may compare it with that produced by a real phase 
function, and adjust a and accordingly. Then we may use the results outlined 
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above to estimate the multiply - scattered contribution from such a phase 
function. Comparisons with second and higher order contributions are also 
possible. This should increase our confidence in the worth of results obtained 
from a model phase function. 


Amush - Stotts (A-S) Type Approximate Method 

In order to extract analytic answers » Arnush (Ref. 38) and Stotts (Ref. 61) 

have expanded P to second order before performing the integration to obtain 0. 

(Series expansion of Q would yield the same result.) Arnush has used Bravo- 

ZhivotovsJtiy's (Ref. 65) sea water phf»se function, Eq. (92a), while Stotts orininall 

used a Gaussian phase function, Eq. (83) and more recently the sea water phase 

function. This approximation is sufficient to provide the correct values for 

2 

both P(z,®), and <r >. 


We start by re-writing the definition of as follows: 


n - 2TI 0) oz (nz) 
o o 


“ 2TI 0 ) oz (nz) 

o 


■t 

■rr 


P(t) dt 


j (ti/;) P(ip) di|» dt 


(104) 


Expanding as a power series leads to 

n = (D oz (1 - n^z^ ^^/12 ...) (105) 

o o 

T 

where is defined by Eq. (83a) . 

It is complicated, but reasonably straightforward to obtain the 
following expression for the beam spread parameter: 

<r^> ^ oz^ + z^/0^ + (106) 

«5 O 

Ignoring the higher order terms in Eq. (105), we may insert this 
expression into Eq. (77') to obtain 
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(107) 


N(z. r) - exp [- (1 - t 



/ n <r^> 


Similarly, integration of Eq. (78) leads to 


P( 2 , R) « F e 
o 

Taking t)ie limits 6, Y 

F(z, R) « F e 
o 

T 

i.e. 0) = e ° 


(1 - U) )T r 2^^ 2n 

° ^ J 


■* “» , we may re-express Eq. (108) in 


(1 - OJ )T 
o 




exp (- 3G /T^) 
s 


•] 




exp {- 3 G / T 




(108) 

more familiar terms 

(108/) 

(109) 


Expansion of to second order in y is equivalent to an asymptotic 
expansion to second order in G or R Thus we may expect this approximation 
to be accurate for large values of G or R. However, its behaviour for small 
values of these parameters is quite different from that of the exact results 
quoted in previous sections. Thus we cannot expect this approximation to 
prove particularly useful for oar problem, as shown later by numerical 
comparisons. In fact, one finds Vbi.ues of A which are less than unity! 

For our problem, of course, we are concerned with small values of R and G, 

and hence we are interested in the behavior of Q for large values of y, i.e., 

o 

its asvmptotic expansion. The phase function model of Eq. (92a) does not have 
an asymptotic expansion, due to the fact that P(0) = «. For the other cases 
we may easily show that 
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(110) 


Vy 

where q - 2ir [ P((|») d4»/a (111) 

Thu« for a Gaussian (Eq. (83)), q - /tI; Eq. (90a) gives q - Ij and Eq. ( 94 a) 
gives q - y B(|^, p + |-) , where B is the beta function. (For W - | example, 

q - 2.) 

For large values of y, the Arnush-Stotts approximation to goes to 
(minus) infinity, and so we cannot expect this approximation to accurately 
predict the power received by a small detector. 

Dolin - Fante (D-F) Approximate method 

Dolin (Ref. 60) and Fante (Ref. 57) have argued that the angular shape of 
the scattered intensity should be a much more slowly varying function than 
P(l/;)« and have thus extracted it from the integral on the right side of Eq. (69) 
After separating the scattered intensity from the unscattered, and Fourier 
transforming, they arrive at the following expressions 

''U . “OZ ^ 

1 (z# 0) - e I (n, Z n) (112a) 

m0 O ^ 


a8 


and 1 (z, n» 0; = I (n 
- ~ o 


n» z n) 

~ it 


dz’ exp [-az'- f X(t, n, Z n) dt) 

'z' 

• P [n (z - z’)] (112b) 
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wh«r« 


(113a) 


19 9 

X(t, n« z n) T «^o V I* n- tn| ♦ (i - u> )o 

tz — • 

X(t,n#e )dt - o 4- (1 - (I) ) o (c • s') (113b) 

jj5» 12 O O 

2 

whera <{« is defined in Eq. (83a) . 

Equations (112) and (113) may now be inserted in Eqs. (77*) or (78) as 
required. Although f is no longer exponentiated, this result is complicated by 
the additional (finite) integration over z*. In the case of a Gaussian phctse 
function, it is possible to reverse the orders of these two integrals, and 
perform that over r|; to give 

fl 2132 

A ■ e - 10 T I dt exp (o) T t - G /(— <u T t + t )1 (114) 

o Jq o 30 

Our calculations show that this approximation is reasonably accurate, 
except in those situations where T is large and G is small. In the Numerical 
Results section, we will compare the predictions of Eq. (114) with those of 
Eq. (88) . As it is not possible to perform any of the integrals for any of 
the other phase function models, we have limited our examination of this 
approximation to the case of the Gaussian phase function. 

NUMERICAL RESULTS 

In this section, we shall present some typical results based on our exact 
formulation and the Arnush-Stotts type approximate method. We shall examine 
4 phase function models: Gaussian (Eq. 83), both exponential models (Eqs. 90a 

and 92a) , and the binomial model with (j - 3/2 (Eq. 94a) . To simplify discussion 
we shall refer to the phase function model of Eq. (90) as the exponential model 
and that of Eq. (92) as the sea-water model. 
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We start by examining the phase functions themselves. It is« of coursc« 
much simpler to plot the normalized phase functions, P, rather than P, where 
P is defined by 

P • 2tt P / (115) 

Unlike P, P is now a function of only one variable, u<)i. In Figure 60, we plot 
P against a<{f, for C ^ a<{* ^ 3.5. We noted earlier that the sea water and binomial 
models have identical values of the parameter q (Eqs. 105). From Fig. 60 we 
see that these two models are quite close over a wide range of values of <3^. 

The next function we examine graphically in (1^. In Figure 61 we plot 
n^/T^ against y for all 4 phase functions, as well as the Arnush-Stotts approxi- 
mation. From this log-log plot, the asymptotic behaviour of the 4 phase functions 
is apparent, particularly that of the sea-water phase function, which has no 
asymptotic expansion. Although it is not obvious from this figure, the curve 
for the binomial phase function actuti.lly lies slightly above the sea-water curve 
for y values less than about 3. Finally we note that for y values greater than 
2, the results obtained by the Arnush-Stotts approximation differ markedly 
from those by our exact formulation, rapidly approaching large negative val^^as 
for y greater than 4. 

We turn now to a discussion of the amplification factor, A, and the power 

received, P, as predicted by these 4 models, and also the Arnush-Stotts type 

approximation. We have evaluated both A and P for G between 0.01 zmd 1.0, and 

T between 0.5 and 15.0. (Hiroughout, we have assumed 0 , Y °°*) 
s 

In the Appendix to this report, we have included a listing of the FORTRAN 
program used to generate this data, along with a brief explanation and sample 
output. 
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Xn Figure 62 we plot A against G for a series of values* for the Gaussian 
model. In Figure 63 we plot P against <assuming unit incident power) for a 

series of values G| again for the Gaussian model. Also shown on this plot is 
the transmission* T* which represents the power that would be received if all 

scattered light was lost. These two graphs clearly indicate the important 
role that forward scattering can play in the detection of transmitted beams* 
especially for optical thicknesses of the order of 10 or higher. 

In Figure 64 we plot A as a function of G for ■ 4.0 and 10.0* for all 
four model phase functions* for our formulation as well as the Arnush Stotts 
approximation. In the latter case* one finds values of A which are less than 
unity. Note that the binomial and sea water curves cross for both values 
(see Pig. 61). For large values of G* we see that there is little to choose 
between the four phase function models. 

In Figure 65, we plot (A - 1.0) against G in log-log form for T ■ 1.0* in 
order to emphasize the linear relationship implied by Eq. (A3) in the Appendix. 

We see that for G less than 0.3* the integral term in Eq. (A3) makes a negligible 
contribution. In Figure 66* we plot (A - 1.0) against G for T ■ 5.0. Here we 
see', that the integral term in Eq. (A3) is starting to make a contribution. Also 
in this graph we have included the Arnush-Stotts and Dolin-Fante approximation 
results. The Dolin-Fante result was not included in Figure 65 as it could not 
be distinguished from the exact result for the Gaussian phase function. 

CONCLUDING REMARKS 

The propagation of a laser beam in an optically dense medium such as a fog* 
dust storm* or smoke is a problem of growing importance* both for communication 
and detection purposes. Although such dense media lead to a significant attenua- 
tion of the primary beam* much of the scattered radiation may still be found 
close to the beam axis, and will thus be available for detection by a suitable 
detector. 
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In thii report, we have examined the spreading of a laser beam using the 
small-angle-scattering approximation to the equation of transfer. This 
approximation appears eminently suited for the study of beam propagation in 
fog, dust, or smoke media, where the scattering phase function is highly 
anisotropic. As well as the standard Gaussian model phase function, three 
other model phase functions have also been examined. The Gaussian functional 
form was used to describe the initial beeun spread and profile, although the 
analysis is somewhat simrler (and the resulting expressions tidier) if the 
limiting case is taken. 

All the numerical results presented in this paper have been based on the 
assumption of a narrow, collimated beam. We may remark, however, that the 
results and expressions presented in this report (e.g., Eq. 86) for the finite 
beam may be applied with full generality. 

We have also examined a number of approximations which have been used to 
further simplify the expressions we have derived. The Arnush-Stotts approxi- 
mation is quite suitable for use in the asymptotic regions, at large distances 
from the beam axis. However, the behavior of the solutions close to the beam 
axis is governed by the pareuneter g, the zeroth moment of the phase function. 
This moment weights the contribution from scattering through very small angles 
far more highly than does the parameter X » scattering angle, or 

third inoment. In fact, for small R (i.e., small G) , one may expand Eq. (88) to 
first order (c.f., Eq.l03) 

A(T ,G) = 1 + T g G + . . . 
s s 


( 116 ) 


Finally, wc may remark that the method proposed by Tam and Zardecki makes 
a useful contribution by providing a connection between realistic (Mie) phase 
functions, and the parameters which must be uoed in the model phase functions 
used in this report. Further work on the applications of this method is 
recommended. 
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NORMALIZED PHASE FUNCTION. 2ir P/a 



a 


FIGURE 60. Normalized phase function P vs. x4* four model 
phase functions. 162 






GEOMETRY FACTOR, 6 


FIGURE 62 . Amplification factor A 
Gaussian phase function. 
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vs. geometry factor G for the 




RECEIVED POWER, P/R 


FIGURE 63. Normalized power received vs. scattering optical thickness 
for the Gaussian phase function. 
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AMPLlFiCATiON FACTOR, A 



FIGURE 64. Amplification factor A vs. geometry factor G for two 
values of Tg/ results are shown for four model phase functions and their 
Arnush-Stotts approximations. 
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P**, iiuht .unpJifio.ition /actor (A - J.O; vs. geometry 

farf»v» (»;* f»M' 1 •• l.O. 


01 - 


100 
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FIGURE 66. Subtracted amplification factor (A - 1.0) vs. geometry 
factor (G) for T ® 5.0. 
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7 . CONCLUDING REMARKS 


Zn this report* we have discussed some new methods, formulations and 
results for the con^lex subject of MS effects on laser beams traversing 
dense aerosols, particularly in connection with topics such as remote 
sensing of time varying aerosol size distributions under MS conditions 
and fast table search methods for retrievals of aerosol size distributions 
undergoing microphysical processes. Conclusions based on work done thus 
far and recommendations for future research are given as follows. 

1. The use of one or two-term analytic models, with each term a 
Modified Gamma Distributions (MGD) or an Inverse Modified Gamma Distribution 
(IMGD) provide reason£d>ly good least-squares fits to the respective unimodal 
and bimodal experimental size distribution data. The selection of the 
analytic model most appropriate for a particular set of experimental size 
distribution depends to a large extent on the asymptotic behavior of data 
plotted on a log-log graph. Details are discussed in Reference 21. 

2. In order to handle a large number of sets of multispectral extinction 
data, to obtain aerosol size distributions, one needs to develop fast retrieval 
techniques. In Sections 4.1, 4.2 and 4.3, we discussed in detail fast table 
search (FTS) methods for two-parameter, three-parameter and three-parameter 
bimodal models, respectively. Comparison of the retrievals obtained by the 
FTS methods and the NLLS method clearly show that the former are 5 to 10 times 
faster than the latter. However, it is strongly recommended that work be 
continued on further optimizing the FTS method with the goal of developing 
automated, on-line size distribution retrieval systems. 

3. Theoretical models have been developed for understanding the 
effects of the dynamical and microphysical processes — sedimentation under 
gravity, coagulation and evaporation/growth — on the temporal behavior of the 
attenuation of visible and IR laser beams traversing dense aerosols. The 
models provide us with some understanding as to the size range in which three 
processes^ — separately and in combinations — dominate the temporal variations 
in aerosol size distribution and, therefore, in aerosol optical extinction. 

Only a few simple cases were studied numerically. To the best of our 
knowledge this is perhaps the most systematic and detailed numerical model- 
ing effort that has been performed on this topic. However, further work is 
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needed and is strongly recommended in order to acquire a deeper understanding 
of the role of these microphysical processes on the attenuation and scattering 
of laser beams traversing the natural environment of hazes, fogs and clouds. 
More work is particularly needed for the case of stirred settling. 

4. In this paper, we developed the formulations for solving the 
radiative transfer equation (RTE) for laser beams traversing scattering 
media in the small-angle approximation and for making the results tractable 
to numerical computations by four different approaches: 

a. Tam and Zardecki exact approach 

b. Arnush-Stotts approach 

c. Dolin-Fante approach 

d. Our exact approach 

The formulations for the four quantities of interest, ncimely, the beam 
radiance, irradiance, detected power and beam spread have been given. 

Numerical computations for the detected power were performed using our 
approach and that of Arnush-Stotts and the two sets of results were compared, 
showing ranges of parameters for which the Arnush-Stotts approach is 
valid. For irradiance computations, our approach differs from that of Teun 
and Zardecki, in that it can numerically treat the case of non-Gaussian 
phase-function. However, the computations have not been performed. It is 
recommended that further work be performed to solve the radiance and 
irradiance of non-Gaussian phase functions by these two approaches. 

5. It should be noted that the small-angle approximation is valid 
for highly forward peaked phase functions. It is, therefore, recommended 
that approaches be developed to deal with the case of MS in laser beams 
traversing small size aerosol particles, with broad phase fun-tions. 

6. It is recommended that the work be done on determining the 
importance of various orders of scattering as a function of the optical 
depth . 

7. It is also recommended that experiments and numerical computations 
be performed to determine the effects of the detector field of view on 
optical extinction measurements. 
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8. It is further recommended that measurementB and theoretical 
computations be performed to inver^tigate the effects of beam diameter on 
optical extinction measurements. 

9. In addition, it is recommended that the above investigation be 
repeated for the case of backscattering, which is of great importance 
for one**cnded electro-optical systems such as Lidars, DIAL and Laser 
Doppler Velocimeters . 
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9 . APPENDIX; COMPUTATIONAL DETAILS 


The nvunerical results presented in Figs. 62 to 66 wore obtained from a 

relatively simple computer program, consisting of less than 100 executable 

statements. A full listing, and partial output, are included in this appendix. 

The task of this program is to evaluate the numerical integrals in 

Eqs. (88), (90d) , (92d) and (95'), for a series of values of between 

0.5 and 15.0, and a series of values of G up to 1.0. For comparison, the 

Arnush-Stotis approximation, Eq. (108), is also computed. We have included 

the full results for T values of 4.0 and 10.0, which may be read in conjunc- 

s 

tion with Fig. 64. 

Although infinite integrals of this type are often handled by Gauss- 
Laguerre quadrature, this method was found wanting, due to the oscillatory 
nature of the integrand. Instead we have employed Simpson's rule, up to a 
finite cut off, allowing for the remainder of the integral by the following 
result: 

If 

f(X) “ 4>(X> for X ” X 

and 

^ 4>(X) dx ■ is known. 


then 



f(X) dX 


^ + f If(X) - <I>(X)1 dx 
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Tu apply this result, we note Eq. (110), and choose X such that 


e ° « 1 ♦ q G/X , X ^ X 


Thus 


r 


Q 

J, (X) e ° dx “ 1 + T q G + 


fX 




Q 

J^(X) (e ° - 1 - q G/X) dx (A3) 


Due to the wide range of values of G which we have used (three orders 
of magnitude) it is necessary to vary the step size accordingly. Thus we 
have used a step size equal to G, up to G ■ 0.22. A step size larger than 
this is unwise, due to the variation in the term. Thus, when we reach 
G * 0.22, a larger set of values is computed and stored, to be used for 
the remaining values of G. 

As pointed out above Eq. (110) , q is undefined for the sea water phase 
function, so we have set q ■ 0 in this case. As a result, we are forced to 
choose a considerably higher value of X in order to satisfy Eq. (A2) . 

In fact, if the sea water phase function was dropped from consideration, 
the time (and cost) of these calculations would be cut at least in half. 

As it is, the results we have obtained for this phase function must be 
considered distinctly lass accurate than the others. 
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